\section{Introduction}
The physics of atomic clusters has gained broad importance as a domain of study of new physical objects, often termed as super-atoms, over the past few decades~\cite{jena2018super}. Such clusters are aggregates of atoms, ranging from two to a few thousand of atoms, and  provide the platform to study a new phase of matter intermediate to atoms or molecules and solids ~\cite{RevModPhys.58.533}. These systems also offer valuable opportunities to understand how the bulk properties emerge from their individual constituent atoms~\cite{stace1988small}. Possibilities of designing new class of materials with tailor made characteristics using clusters as building blocks, instead of atoms, have rendered this field diverse and vibrant~\cite{castleman2009clusters}. The extensive research on clusters over the years has revealed many unusual properties which are of interests to the field of physical and chemical sciences, material sciences, and biological sciences, making these explorations interdisciplinary in nature~\cite{dinh2013introduction}.

Many intriguing phenomena are associated with the photoresponse spectrum of atomic clusters~\cite{kreibig2013optical}. These include phenomena, such as, plasmon resonances~\cite{xia2009,madjet2009,ekardt1985}, Fano type resonances~\cite{miroshnichenko2010fano}, and diffractive modulations~\cite{jaenkaelae2011} in the photoelectron signal due to the largely well-defined cluster edges~\cite{madjet2001}. One primary focus of the present work is the photoionizing response of the plasmon resonances - resonances that form due to the collective oscillations of the valence electron cloud. Excitation of a lower energy giant surface plasmon resonance is known to be a prominent feature in the photoresponse of alkali metal clusters below the ionization threshold. This feature in metal clusters has applications in nano-optical devices, chemical and biological sensing, and bio-medicine etc~\cite{mirin2009,loo2005,liao2006}, besides its eminent role as ``spectral laboratories" to assess many-electron effects. Therefore, developing detailed and accurate understanding of the origin, the underlying dynamics, and related observable effects of this resonance are a matter of significant priority. 

In addition, the presence of a higher-energy volume like plasmon makes the photo-spectrum of atomic clusters~\cite{xia2009} more robust compared to the spectrum of the corresponding bulk metal. It is well known that due to the translational invariance, the coupling of the volume oscillation in bulk to light is not feasible~\cite{ferrell1958predicted}. However, the situation is different for the case of finite systems like metal clusters due to the broken translational invariance that may cause boundary reflections of the plasma wave. As a result, optically active volume plasmon resonance can also become possible. The presence of such a volume plasmon excitation is observed for the first time in the experimental work of Xia et al.~\cite{xia2009}. Their photon depletion measurements of Na$_{20}$ and Na$_{92}$ clusters showed a broad volume plasmon resonance mounted on the decay ridge of giant surface plasmon and peaking slightly above 4 eV. The tail of this combined structure was found extended to the ionization region. The present calculation is motivated to address this spillover part of the plasmon structure closely above the ionization threshold. In studying this, we also predict a universal minimum at the decay end of the structure from an interchannel coupled phase-coherence effect arising from many-electron correlations. Furthermore, an analysis of the many-body induced self-consistent field potential at the plasmon spillover energies point to a correlation driven attractive force that suggests possible time delays of the emerging photoelectron. 

The dynamical response of the Na clusters to the external electromagnetic radiation is calculated using a linear response density functional theory (DFT) approach. A fairly competent method belonging to this class is the jellium based time-dependent local density approximation (TDLDA)~\cite{onida2002electronic}. In the past, application of TDLDA to atomic clusters is found to be successful in explaining the collective phenomena occurring at low energies~\cite{madjet2009} and also the diffractive oscillations present at high energies~\cite{madjet2001}. The ease and transparency of a jellium based model enable physicists to interpret the key physics that determines the dynamics. A recent study on fullerenes indicated that a better agreement with experiment can be obtained by using one such DFT approach with Leeuwen-Baerends (LB94) exchange-correlation functional that produces the correct ground state asymptotic properties~\cite{choi2017effects}. Another study has tested the efficacy of the method for Na clusters~\cite{shaik2020}. Hence, this method is chosen in the present work.
  
The paper is organized in the following way. A basic description of the theoretical methodology is provided in Section II which has two parts. The details of ground state structure in the spherical jellium formalism along with a brief account of the LB94 parametrization scheme is provided in Sec. II A. A brief description of the method that incorporates electron correlations in response to the radiation is given in Sec. II B . Section III constitutes the results and discussions of: photoionization cross sections and its comparison with experiment (Sec. III A), the photoelectron oscillator strength (Sec. IIIB), the many-body correlation induced potential in the spillover region of plasmon (Sec. III C), and the prediction of a ``correlation minimum" at the high-energy end of the plasmon structure (Sec. III D). Finally, Section IV concludes the study.

\section{Theoretical Methodology}

\subsection{Ground states of Na$_{20}$ and Na$_{92}$}

To investigate the ground state electronic structure of Na$_{n}$($n$=20 and 92), a DFT approach is adopted. For such systems of closed-shell configurations the jellium model serves as a very good approximation. In this model, the jellium potential ($V_{jel}(\bm r)$) replaces the ionic core of 20 and 92 Na$^{10+}$ ions, respectively for Na$_{20}$ and Na$_{92}$, by potentials constructed after homogeneously smearing their positive charges into jellium spheres. The radius of each cluster is determined by the number $n$ of ions present in the system. The radius of Na$_{20}$ is calculated to be 10.67 a.u.\ and that of Na$_{92}$ to be 17.74 a.u.\ using the formula  $R_{c}$=$r_{s}n^{1/3}$, where $r_{s}$ is the Wigner-Seitz radius (3.93 a.u.) of a Na atom. The Kohn-Sham equations for $n$ delocalized valence electrons, a 3s$^1$ electron from each Na atom, are solved to obtain the ground state structures of Na$_{20}$ and Na$_{92}$. It is to be noted that to match the valence ionization thresholds with the experimental values~\cite{Chandezon}, suitable constant pseudo potentials are added.  
\begin{figure}[ht]
\includegraphics[width=\linewidth]{ground.eps}
\caption{a) Ground-state radial wavefunctions for HOMO and HOMO-1 levels calculated for Na$_{20}$ and Na$_{92}$. b) Ground-state radial potentials and energy levels calculated for Na$_{20}$ and Na$_{92}$.}
\label{gpot}
\end{figure}

In terms of the single-particle density $\rho(\bm r)$, the ground state self-consistent field LDA potential reads as,
\begin{equation}\label{lda-pot}
V_{\scriptsize \mbox{LDA}}(\mathbf{r}) = V_{\mbox{\scriptsize jel}}(\mathbf{r}) + \int d\mathbf{r}'\frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} + V_{\scriptsize \mbox{xc}}[\rho(\mathbf{r})],
\end{equation}
where the second and third terms on the RHS are respectively the direct and exchange-correlation (xc) potentials. Approximate form for $V_{\scriptsize \mbox{xc}}$ has to be used since its exact form is unknown. In the present work, we approximate the $V_{\scriptsize \mbox{xc}}$ with the LB94 functional which provides an accurate asymptotic description of the ground-state properties. This functional belongs to a gradient corrected class of approximation by van Leeuwen and Baerends, parameterized in terms of the reduced density and its gradient $\nabla\rho(\bm r)$~\cite{van1994}, to be given by,\\
\begin{equation}\label{lb94}
V_{\mbox{\scriptsize LB}} = -\beta{[\rho({\bf r})]^{1/3}}\frac{(\xi X)^2}{1+3\beta\xi X \sinh^{-1}(\xi X)},
\end{equation}
where $\beta$=0.01 is a fitting parameter and $X = [\nabla\rho]/\rho^{4/3}$. The parameter $\xi$ arises due to the change of system from spin-polarized to spin-unpolarized form ~\cite{oliver79}. This scheme of gradient correction of density is derived from a gradient expansion series that naturally eliminates self-interactions and produces the correct asymptotic behavior.

\begin{table}[ht]
\caption{Binding energies (BE) of HOMO and HOMO-1 levels of Na$_{20}$ and Na$_{92}$ in the harmonic oscillator notation.}
\label{tab1}      
\begin{tabular}{c|c cccc}
\hline\noalign{\smallskip}
    &    & BE$_{\mbox{\scriptsize HOMO}}$ (eV) &    & {\smallskip} BE$_{\mbox{\scriptsize HOMO-1}}$ (eV)\\
\noalign{\smallskip}\hline\noalign{\smallskip}
Na$_{20}$  &  & $2s$ {\smallskip} $-3.75$ & &$1d$ {\smallskip} $-4.20$      \\
Na$_{92}$ &     & $1h$ {\smallskip}  $-3.47$ & & $3s$  {\smallskip} $-3.57$  \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

The calculated binding energies of the highest occupied molecular orbitals (HOMO) and (HOMO-1) of these systems,  of 2s and 1d character respectively for Na$_{20}$ and 1h and 3s respectively for Na$_{92}$, are given in Table~\ref{tab1} and their radial wavefunctions are shown in Fig.~\ref{gpot}(a). It can be seen that HOMO level of  Na$_{20}$ has a single node whereas Na$_{92}$ is nodeless. The situation reverses for HOMO-1. As a result, the radial structures of HOMO and HOMO-1 wavefunctions significantly differ in these two systems. Fig.~\ref{gpot}(b) shows the ground-state radial potentials of Na$_{20}$ and Na$_{92}$ clusters and the level energies.
 
\subsection{Dynamical Response}
The dynamical response of the cluster subjected to the external dipole field `$z$' can be calculated by employing a linear-response DFT scheme, known as time-dependent LDA (TDLDA) modified by LB94~\cite{madjet2008photoionization}. The dipole field can induce a frequency-dependent complex change in the electron density $\delta \rho$~\cite{petersilka1996}, which, within a linear-response perturbative framework, can be written as:
\begin{equation}\label{cor-dens}
\delta \rho (\mathbf{r}^{\prime}; \omega) = \int \chi (\mathbf{r}, \mathbf{r}^{\prime}; \omega)
z  d\mathbf{r},
\end{equation}
where $\chi$ is the full susceptibility of the system which incorporates all the dynamical electron correlations. The above equation can be re-written in terms of independent particle (IP) susceptibility, $\chi_0$, and the complex field $\delta V$ as, 
\begin{equation}\label{cor-dens-1}
\delta \rho (\mathbf{r}^{\prime}; \omega) = \int \chi_0 (\mathbf{r}, \mathbf{r}^{\prime}; \omega)
\delta V(\mathbf{r}; \omega)  d\mathbf{r},
\end{equation}
where $\chi_0$ is constructed by the ground-state single-electron orbitals $\phi_{nl}$ and energies $\epsilon_{nl}$~\cite{feibelman1975microscopic}. The IP susceptibility $\chi_0$ is related to $\chi$ by the matrix equation: 
\begin{eqnarray}\label{chi-chi0}
\chi = \chi_0[1-(\partial V/\partial \rho)\chi_{0}]^{-1}.
\end{eqnarray}
The total field, $\delta V$, is given by the following relation: 
\begin{equation}\label{comp-pot}
\delta V(\mathbf{r}; \omega)=z+V_{ind}(\mathbf{r}; \omega),
\end{equation}
where
\begin{equation}\label{V_ind}
V_{ind}(\mathbf{r}; \omega) = \int \frac{\delta \rho (\mathbf{r}^{\prime}; 
\omega)}{|\mathbf{r}-\mathbf{r}^{\prime}|}d\mathbf{r}^{\prime}+\Big[\frac{\partial V_{XC}}{\partial \rho}\Big]_{\rho = \rho_{0}} \delta \rho(\mathbf{r}; \omega).
\end{equation}

Using the matrix inversion method ~\cite{BERTSCH1990}, Eq.~\ref{chi-chi0} is finally solved for $\chi$ which in turn is used for obtaining $\delta \rho$ and thereby $\delta V$ by making use of Eq.~\ref{cor-dens} and Eq.~\ref{comp-pot} in a self-consistent way.


In this formalism, the photoionization cross section corresponding to a bound-to-continuum dipole transition $n\ell\rightarrow k\ell^\prime$ is then calculated using
\begin{equation}\label{cross-pi}
\sigma_{n\ell\rightarrow k\ell'} \sim |\langle k\ell'|\mbox{$\delta V$}|n\ell\rangle|^2.
\end{equation}
It is clear from the above equation that in addition to the external perturbation $z$, the calculation involves the complex induced field $V_{ind}$ driven by the many-electron correlations. Obviously, setting $\delta V=z$ yields the independent particle LDA cross section  that ignores correlations. A comparison of LDA and TDLDA, therefore, easily possible in this method to study the role of many-electron effects in the photoionization process. 
\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig1.eps}
\caption{TDLDA and LDA cross sections of Na$_{20}$ (a) and Na$_{92}$ (b).}
\label{lda-tdlda}
\end{figure}

\section{Results and Discussion}

\subsection{Photoionization cross-sections and comparison with measurements}

Total photoionization cross-section of Na$_{20}$ and Na$_{92}$ clusters, calculated in TDLDA+LB94 and the corresponding LDA results, are shown in Fig.~\ref{lda-tdlda}(a) and Fig.~\ref{lda-tdlda}(b) respectively. As seen, the TDLDA cross section profiles, as a function of photon energies, strongly differ from the LDA profiles by (i) significant enhancements near the ionization threshold region, (ii) the presence of autoionization resonances, and (iii) the appearance of a minimum at about 7.20 eV for Na$_{20}$ and 6.22 eV for Na$_{92}$. The near-threshold enhancement in TDLDA cross sections combines the remnant of the giant surface plasmon resonance, that flourishes below the ionization threshold, and a major portion of the broader volume plasmon resonance, both emerging from the collective electronic motions when subjected to an external electromagnetic field. %
This resonance spillover from the discrete to continuum spectrum is embedded by a host of narrow spikes, which are the autoionization resonances resulting from the degeneracy of ionization channel with the inner-level single-electron discrete excitations~\cite{fano1961effects}. As seen, these features are completely missing in the LDA predictions which neglect the electron correlations. Structures in the LDA profiles in this near-threshold region are due to the gradual openings of inner-level ionization channels.  
\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig2.eps}
\caption{TDLDA cross sections along with experimental data from ~\cite{xia2009} for Na$_{20}$ (a) and Na$_{92}$ (b). Lorentzian fits to the experimental cross sections are also shown.}
\label{tdlda-expt}
\end{figure}

The existence of the volume plasmon resonance and the extended tail of the surface resonance to the ionization part is experimentally seen in the work of~\cite{xia2009}. Comparisons of our calculations with these measured data is shown in Fig.~\ref{tdlda-expt}(a) for Na$_{20}$ and Fig.~\ref{tdlda-expt}(b) for Na$_{92}$. The data covers the range of 2 eV to 5.64 eV. Following Ref.\ \cite{xia2009}, fits to the experimental data as a sum of four Lorentzian profiles for Na$_{20}$ and a sum of three for Na$_{92}$ are also shown in Fig.~\ref{tdlda-expt}. The lower energy Lorentzians in both cases represent the giant portion of the plasmon below the first ionization threshold that our photoionization results can not access. The present TDLDA calculations span the ionization part of the spectrum only and, as noted, overlap very well with the higher energy Lorentizian fits used for both the clusters. The peak energy position (E$_0$) and width ($\Gamma$) of the fourth Lorentizian curve for the higher energy plasmon resonance (HEP) for Na$_{20}$ and the corresponding position and width for the third Lorentzian curve for Na$_{92}$ are shown in Table~\ref{tab2}. It is clear from the figure that the steady background parts of TDLDA cross sections agree well with the experiment describing the spillover of the plasmon resonances. As the size of the cluster increases the peak positions of the volume plasmon shifts to higher energies (blue shifted). Curiously, this size effect is found to be opposite in the case of fullerenes (C$_{60}$ and C$_{240}$) where peak positions shifted to lower energies as the size increased~\cite{choi2017effects}. It is to be noted that the narrow autoionization resonances are missing in the experimental data. This is due to the finite temperature effects of the metal clusters in experimental conditions which lead to the coupling of electronic motion with the temperature induced vibrational and rotational modes~\cite{bertsch1989thermal}, as was also found and discussed earlier for C$_{60}$~\cite{madjet2008photoionization}. 

\begin{table}
\caption{Resonance positions ($E_o$), FWHM ($\Gamma$) of the higher energy plasmon (HEP)~\cite{xia2009}.}
\label{tab2}      
\begin{tabular}{c  c  c}
\hline\noalign{\smallskip}
              & $E_o$ (eV)   & $\Gamma$ (eV) \\
\noalign{\smallskip}\hline\noalign{\smallskip}

Na$_{20}$ HEP  & 4.04  & 1.19     \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Na$_{92}$ HEP & 4.20  & 1.16   \\   
\noalign{\smallskip}\hline
\end{tabular}
\end{table}



\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig3.eps}
\caption{Accumulative oscillator strengths for Na$_{20}$ (a) and Na$_{92}$ (b) where the theoretical results are calculated using the present TDLDA+LB94 cross sections and the experimental OS are calculated using fitted values to the data in~\cite{xia2009}.}
\label{OS}
\end{figure}

\subsection{Oscillator strength} 

From the experimental data shown in Fig.~\ref{tdlda-expt}, it is quite clear that the bulk of the absorption oscillator strength (OS) is consumed in the excitation part of the spectrum below the ionization threshold. In order to get a quantitative measure of the fraction of the OS exhausted in the ionization process, we introduce an ``accumulative" OS as a function of photon energy given by the following relation:   
\begin{equation}\label{OSt}
OS_{t}(E_{f})=OS_{b}+\int_{E_{th}}^{E_{f}} \sigma(E)dE, 
\end{equation}
where OS$_{b}$ is defined as the baseline OS that is calculated by integrating the experimental cross section~\cite{xia2009} from the starting photon energy of the measurement to the theoretical ionization thresholds (E$_{th}$). OS$_{b}$ therefore corresponds to the oscillator strength exhausted by the plasmon resonances below the threshold. The second term in Eq.~\ref{OSt} corresponds to the incremental addition to the OS$_{b}$, cumulatively with energy, due to the plasmon spill to the ionization region from E$_{th}$ to E$_{f}$. Here E$_{f}$ will vary from E$_{th}$ onwards allowing us to calculate OS$_{t}$ as a function of E$_{f}$. Note that until 5.64 eV, the photon energy corresponding to the last experimental datum, E$_{f}$'s are the various energies for which experimental measurements are available. In Eq.~\ref{OSt}, $\sigma(E)$ corresponds to TDLDA total cross section for the calculations of theoretical OS$_{t}$. For calculating the experimental OS$_{t}$, on the other hand, fitted cross section data between two successive measurements are used.   

Figure~\ref{OS} shows OS on a zero-to-one scale, because it represents OS per atom, that is the total OS divided by $n$. For Na$_{20}$, the OS$_{b}$=0.5230 and for Na$_{92}$, OS$_{b}$=0.4352. These suggest that about 52\% of OS is exhausted below the ionization threshold of Na$_{20}$, while it is about 44\% for Na$_{92}$. The implication is that more fraction of electrons are available to participate in the ionization process for the larger cluster. Hence, the size of plasmon spill increases for  Na$_{92}$. The rise of the theoretical OS curve above the experimental curve is due primarily to the narrow autoionizing resonances in the theoretical spectrum which add strengths. This accounts roughly about 20\% of OS for Na$_{20}$, while it is 30\% for Na$_{92}$. The experimental measurements cannot access such narrow single particle resonances due to the finite temperature effects~\cite{bertsch1989thermal}; similar distinction was also noted for the plasmon resonances in C$_{60}$~\cite{madjet2008photoionization}. However, the missing part of the oscillator strengths in the measurements may be distributed over higher photon energies, in which case it is expected that the two curves may eventually intersect at higher energies. One comment should be made about the measurements though. The energy deposited into the cluster will not only dissipate through ionization that we calculate, but also be shared with thermalization and evaporation channels as well. Even though the corresponding branching ratio is unknown, the ionization is known to be a much faster process (in attoseconds) compared to thermalization/evaporation (in femtoseconds). Therefore, one may reasonably assume that the main contribution to the measured data above the threshold is likely from the ionization process. 
\begin{figure*}[!ht]     
\includegraphics[width=16cm]{scf.eps}
\caption{The real and imaginary part of radial self-consistent field potential, $V_{ind}(r)$, for Na$_{20}$ (left panels) and Na$_{92}$ (right panels).}
\label{scf}
\end{figure*}

\subsection{Collectivity induced field}

TDLDA dipole matrix element $ \langle k\ell'|z+V_{ind}|n\ell \rangle$ in \eq{cross-pi} requires knowledge of $V_{ind}$, which is the complex induced field driven by electron correlations and is singularly responsible for the plasmonic enhancement in the cross section. The behavior of the real and imaginary parts of $V_{ind}$ across the collective resonance region is well known~\cite{wendin73,zangwill80}. Im($V_{ind}$) characteristically shows a well-type shape across the energy range of the resonance where the minimum of the well occurs near the energy of the resonance peak. On the other hand, Re($V_{ind}$) executes an oscillation by switching the sign over this range where it sluices through the zero at the resonance peak. These two distinct behaviors can be combined in a unified picture. Im($V_{ind}$) has a predominant collective character, while Re($V_{ind}$) represents effects of the external field. Therefore, as the resonance builds with increasing energy and approaches its peak, the effect of external field reduces while the collective motion grows. At the peak, the effect of external field is negligible, where the collective response dominates. At the waning part of the resonance this trend reverses.

In the current study, however, we access the remnant of the plasmon structure spilling over to the photoionization channel which is fragmented by many single electron narrower resonances mixing coherently with the plasmon effect. This coherence indicates that there are interferences between the single electron Auger and collective plasmon. So the simple mechanism expressed above becomes complex, the induced field becomes structured, but an interesting general trend should remains that we now explain. 

Figure~\ref{scf} shows the 3D plots of real and imaginary part of radial induced field, $V_{ind}(r; \omega)$, for Na$_{20}$ and Na$_{92}$. All the plots show many  smaller structures, which are because of single-electron resonances that become relatively stronger as the plasmon spill gradually weakens with increasing energy.  However, there appears wider well-shapes for Im ($V_{ind}$) of both systems. This broader range of negative fields suggests the emergence of an attractive field as a result of the electrons' collective dynamics in these plasmonic energies. Consequently, the liberating photoelectrons will experience a resistance against their exit and will likely slow down. This can lead to observable time-delayed emissions of electrons~\cite{schultze} in the plasmon spillover range. Particularly, this may find relevance in the context of the electron's intrinsic Wigner-type delay properties~\cite{magrakvelidze2015attosecond}. Time delay studies in photoemission belong to a contemporary field of interest based on methods of attosecond photoemission measurements in RABITT or streaking spectrometry~\cite{Pazourek}. Similar attractive force, driven by electronic collective interactions, has earlier been predicted at plasmon photoionization energies of C$_{60}$~\cite{magrakvelidze2020}. 
\begin{figure*}[ht]
\includegraphics[width=\linewidth]{fig4.eps}
\caption{Subshell cross-sections for Na$_{20}$, (a) and (c), and Na$_{92}$, (b) and (d), compared with corresponding total cross-section obtained using LDA (left-panels) and TDLDA (right panels).}
\label{subcross}
\end{figure*}

\subsection{Correlation minimum}

As seen in Fig.~\ref{lda-tdlda}, the total photoionization cross sections undergo oscillations both in LDA and TDLDA starting from energies where the plasmon structure begins to fizzle. The oscillations far above the plasmon region are well understood and are attributed to the diffraction of the photoelectron waves from the edges of the cluster~\cite{madjet2001}. Since the feature is associated with the cluster geometry, it is seen even in the independent particle LDA results. The merging of the oscillations at higher energies in TDLDA and LDA reflects on the fact that the oscillations have nothing to do with the electron correlations. However, the TDLDA cross sections reveal a rather strong additional minimum occurring at 7.20 eV for Na$_{20}$ and 6.22 eV for Na$_{92}$ which are clearly missing in their LDA counterparts. Incidentally, this feature was also revealed in earlier calculations~\cite{ekardt1985}, but was never properly interpreted. If it were a diffractive feature, it would be present in the LDA results as well. Obviously, that is not being the case suggests that this lower energy minimum must be induced by an electron correlation effect. 

Further insights into this minimum can be harnessed by looking at the individual subshell cross sections. In Fig.~\ref{subcross} are shown the cross sections from the three outermost subshells of Na$_{20}$ and Na$_{92}$ along with the total cross sections obtained from LDA (Fig.~\ref{subcross}(a) and Fig.~\ref{subcross}(b)) and TDLDA(Fig.~\ref{subcross}(c) and Fig.~\ref{subcross}(d)). It is clear from the TDLDA results that in each of the subshell cross sections for a cluster this particular minimum in consideration occurs dramatically at the same low photon energy to create a rather emphatic minimum in the total cross section. But this feature is totally absent in the LDA curves. 

For a subshell, while the creation of a minimum derives directly from the ionization amplitude, its position in energy traces in to the ionization phase. The positions of the minima at high energies are seen to be very sensitive to the ionizing subshells. It was shown earlier that the energy-positions of high energy minima are a function of the phase of the individual subshell LDA transition amplitudes~\cite{frank97}, since the contribution of $\delta V$ [\eq{comp-pot}] is virtually zero at these energies. In addition, at very high energies the contributions to the scattering-phase from even the short-range and Coulomb potentials are negligible. As a result, the final state wavefunction is approximately of the form $\psi_f \approx \cos(kr-\frac{l^{\prime} \pi}{2})$, where $l^\prime = l \pm 1$ being the dipole selected final state angular momentum, `$l$' is the initial state and `$k$' is the photoelectron momentum. Consequently, the squared transition amplitudes differ by a phase of $l^\prime \pi$. This explains why the oscillations in the cross sections of two subshells of angular momenta differing by an odd integer, such as, s and p, p and d, s and h etc.\ are roughly out-of-phase, while those differing by an even integer, such as, s and d, p and f etc.\ are roughly in-phase~\cite{mccune2008unique}. This can be seen in both LDA and TDLDA subshell results at very high energies. At not-so-high energies, these patterns are not seen to be exactly followed due to non-negligible short-range and Coulomb phases. 


However, the energy-concurrent minima in the TDLDA profiles occur at a low enough energy. Even though the diffractive mechanism forming oscillations is emerging, this range is still in the waning part of the plasmon structure. Hence, the contribution from the correlation phase, that is the phase of the amplitude with significant strength of $\delta V$, still dominates over the $l^\prime \pi$ in this range. The correlation phase, being collective in nature, is subshell independent, leading to all subshells experiencing effectively the same phase~\cite{magrakvelidze2015attosecond}. This ensures that the minimum for all the subshells to appear coherently at the same energy. Since this coherence originates from a direct influence of the collective dynamics, we call this feature a correlation minimum, which can, in principle, be accessible by the photoelectron spectroscopy (PES). 

\section{Conclusions}      
In summary, the present jellium based linear response DFT calculation describes the remnant of the giant surface and the bulk of the volume plasmon in the photoionization cross section of Na$_{n}$ ($n$=20 and 92) clusters.  Results show that an appreciable amount of plasmon spillover into the ionization continuum occurs in these systems. The steady background part of the cross section exhibits reasonably good agreement with previous experimental results of the absorption of these clusters. However, a detailed scrutiny of the  oscillator strength calculated from the cross section and its comparison with the strength extracted from the measurements enables us to quantify the contribution of the single electron Auger-type resonances that the theory predicts. Such, rather narrow, resonances shown in the theory, which does not incorporate the temperature and vibro-rotational effects, are missing in the measurements that include these effects. Furthermore, a deeper scrutiny of the many-body induced self-consistent DFT field reveals the presence of an attractive force in the plasmon spillover energy region. We speculate that this force can cause an observable delay in the emission of photoelectrons at these energies, enabling these clusters as interesting candidates for time-delay measurements. The current study further uncovers the presence of a correlation minimum in the cross sections which appears at the waning range of the plasmon resonance structure and which may also attract experimental effort. To this end, the collective-motion spectral properties divulged in the current DFT study of Na clusters are likely quite general and may apply at various degrees of prominence in other metallic clusters. 
\begin{acknowledgments}
The research is supported by the SERB, India, Grant No.\ EMR/2016/002695 (HRV) and by the US National Science Foundation Grant No.\ PHY-1806206 (HSC).
\end{acknowledgments}

