Density Functional Treatment of Photoionization of Sodium Clusters: Effects of Cluster Size and Exchange–Correlation Framework

Abstract

The ground state and photoionization properties of Na${}_x$ (x = 20, 40, and 92) clusters are investigated using a method based on density functional theory (DFT) in a spherical jellium frame. Two different exchange–correlation treatments with the Gunnarsson–Lundqvist parametrization are used: (i) the electron self-interaction correction (SIC) scheme and (ii) the van Leeuwen–Baerends (LB94) scheme based on the gradient of the electron density. The shapes of the mean-field potentials and bound state properties, obtained in the two schemes, qualitatively agree, but differ in the details. The effect of the schemes on the photoionization dynamics, calculated in linear response time-dependent DFT is compared, in which the broader features are found to be universal. The general similarity of the results in SIC and LB94 demonstrates the reliability of DFT treatments. The study further elucidates the evolution of the ground state and ionization description as a function of the cluster size.

1. Introduction

Atomic and molecular cluster physics has emerged as a distinct area of research over the last decades. The area has evolved into an important field by bridging the gap between the atomic/molecular and the bulk domain. This opened new pathways to characterize and control parts of the nano-world. The research that emerged has intertwined physics, chemistry, astronomy, and biology, thus making cluster studies an interdisciplinary topic [1,2,3]. Theoretical models for the description of atomic clusters based on ab initio principles along with the experimental studies have unravelled intriguing features, especially resulting from the interaction of a cluster with light [4,5]. Such photo-induced processes include, for instance, the plasmon resonances [6,7,8,9], Auger-type Fano resonances [10,11], inter-Coulombic decay (ICD) resonances [12,13,14,15], and modulations in the photoelectron intensity due to the diffraction from cluster edges [16,17,18]. In addition to their role as “spectral laboratories” to probe many-electron effects, cluster studies provide impetus for a wide variety of applications, such as, in nano-optical devices [19], in solar energy harvesting [20], and in chemical and biological sensors [21,22,23]. Therefore, testing the efficacy of theoretical models, taking into account the many-body and quantum phenomena, is particularly valuable.

Experiments suggest that the details of the ionic core configuration of metal clusters, such as sodium clusters, play a less significant role in extracting structural and dynamical information [24,25]. In addition, the loosely bound valence electrons in the clusters can be approximated as delocalized and confined within a broad potential well. The jellium model, which makes use of these facts, replaces the metallic ion core by a uniform charge distribution, which provides an electrostatic attraction to the valence electron cloud. The electronic structure may then be determined by applying a mean-field approximation to the interacting electrons that includes static exchange and correlation ($XC$) effects in addition to the direct electron–electron Coulomb repulsion.

An accurate description of the $XC$ potential is crucial in the above model. The Kohn–Sham (KS) density functional theory (DFT) method is known to have some limitations in handling the electron exchange, which is intrinsically a fully non-local attribute, by treating it in a local frame [26]. As a result, the exchange functional cannot perfectly cancel the self-interaction present in the direct (Hartree) term. This causes the $XC$ potential for a finite system to approach zero exponentially as $r\to \infty$ rather than to residual $-1/r$, which is the correct long-distance behavior for a neutral system. Certain structural and dynamical properties are affected due to this inaccurate description of the asymptotic behavior. It leads to errors in the determination of properties, such as the ionization potential and cross-section. In order to overcome this limitation, for instance, one of the following approximate methods can be used: (i) the self-interaction correction (SIC) or (ii) the generalized gradient approximation (GGA).

In method (i), introduced by Perdew and Zunger [27], the self-interactions are subtracted from the potential in the KS equations and iterated until a self-consistent solution is obtained. It improves the asymptotic behavior of the DFT potential, although the resulting KS Hamiltonian becomes state-dependent. This approach was found to be useful for a wide range of compact atomic or molecular systems [28], and especially so for explaining the absorption spectra of alkali metal clusters [29].

Method (ii), which is more intrinsic to the formalism, is within the GGA class and was developed by van Leeuwen and Baerends [30]. In this approach, known as LB94, the issue with self-interaction is addressed by introducing a term that is dependent on the gradient of electron density by using the Becke GGA construction for the modeling of the $XC$ potential. LB94 produces a state-independent potential, thus offering a relatively easy and inexpensive implementation in the computer code. The study of fullerene molecules using LB94 is found to show a somewhat better agreement with the experimental results [31]. However, in a recent ICD study of fullerene plasmon resonance in Na${}_{20}$@C${}_{240}$, the SIC ground state structure complied better with quantum chemical results [15].

The photoresponse of alkali metal clusters energetically below the ionization threshold is dominated by the giant surface plasmon resonance excitation. The photospectra of clusters are more robust than their bulk counterpart due to the existence of a higher-energy volume plasmon [6]. Our previous theoretical study on Na-clusters [32] predicted the spillover of a smaller remnant of the giant surface plasmon and most of the volume plasmon out to the ionization region and showed a reasonable agreement with measurements. It further predicted a feature called the correlation minimum in the ionization spectra of Na clusters. However, in [32], the LB94 scheme was employed while modeling the $XC$ interaction and a linear response time-dependent DFT (LR-TDDFT) formalism was used to calculate the response properties. Since the photoionization spectra are expected to be somewhat sensitive to the $XC$ model, in this study, we compare the effect of the two schemes, i.e., SIC and LB94, on the details of the photoresponse behavior using LR-TDDFT. Furthermore, in order to explicate the evolution trend of the results on the cluster size, three cluster systems are selected for the study.

2. Theoretical Methodology

2.1. DFT Exchange–Correlation Functionals

The details of the method are in line with the framework discussed in Ref. [33]. We adopted a jellium-based DFT approach to explore the ground state electronic structure of Na${}_x$ (x = 20, 40, and 92) clusters in a spherical model. The jellium potential, $V_{\mathrm{jel}}$, replaces the ionic core of N (N = 20, 40, and 92) Na${}^{+}$ ions with a potential created by homogeneously smearing their total charges over a sphere. The radial component of the spherically symmetric potential generated by this distribution is the following:

$$V_{\mathrm{jel}}\left(r\right)=\left\{\begin{array}{c}-\frac{N}{2R_c}(3-{\left(\frac{r}{R_c}\right)}^2),r\le R_c\hfill \\ -\frac{N}{r},r>R_c\hfill \end{array}\right.$$

The radius of each cluster is determined by $R_c$ = $r_s{N}^{1/3}$, where $r_s$ = 3.93 a.u. is the Wigner–Seitz radius of the Na atom. The KS equations for N number of delocalized valence electrons, i.e., the 3s${}^1$ electron from each Na atom, are solved to obtain the electronic structure. It is to be noted that a constant pseudo-potential is added to match the first ionization threshold with the experimentally known values [34].

The ground state self-consistent field DFT potential can be written in terms of the single-particle density $\rho \left(\mathbf{r}\right)$ as,

$$V_{\mathrm{DFT}}\left(\mathbf{r}\right)=V_{\mathrm{jel}}\left(\mathbf{r}\right)+\int d{\mathbf{r}}^{\prime }\frac{\rho \left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+V_{XC}\left[\rho \left(\mathbf{r}\right)\right],$$

(1)

where the second and third terms on the RHS are known as the direct and the $XC$ components, respectively. We consider the following formula to initially parameterize V${}_{XC}$ by using $\rho \left(\mathbf{r}\right)$ [35]:

$$V_{XC}\left[\rho \left(\mathbf{r}\right)\right]=-{\left(\frac{3\rho \left(\mathbf{r}\right)}{\pi }\right)}^{1/3}-0.0333log\left[1+11.4{\left(\frac{4\pi \rho \left(\mathbf{r}\right)}{3}\right)}^{1/3}\right].$$

(2)

Equation (2) qualifies the method as a local density approximation (LDA) approach. By the standard variational technique, it is possible to exactly derive the first term on the right hand side of Equation (2) from the Hartree–Fock (HF) exchange energy of a homogeneous electron system that has a uniform positively charged background. The second term is the so-called correlation potential, which is not accounted for in the HF formalism. As mentioned earlier, the localization of the potential leads to the non-cancellation of self-interactions. A corrective scheme is therefore adopted from the outset to artificially eliminate unphysical self-interactions for each $i\mathrm{th}$ occupied subshell. This leaves the LDA potential orbital-specific, but it approximately captures the asymptotic properties of the electron. We describe this model, referred to as DFT-SIC, in the equation below:

$$V_{\mathrm{DFT}-\mathrm{SIC}}^i\left(\mathbf{r}\right)=V_{\mathrm{jel}}\left(\mathbf{r}\right)+\int d{\mathbf{r}}^{\prime }\frac{\rho \left({\mathbf{r}}^{\prime }\right)-{\rho }_i\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+V_{XC}\left[\rho \left(\mathbf{r}\right)\right]-V_{XC}\left[{\rho }_i\left(\mathbf{r}\right)\right].$$

(3)

DFT-SIC thus mimics two desirable $XC$ functional features: it cancels out the self-interaction part of the Hartree energy and it vanishes for a one-electron system.

An alternative method to correct the $XC$ functional makes use of Equation (2) but refines it further by the addition of a parameterized potential defined in terms of the reduced density and its gradient $\nabla \rho$, as follows:

$$V_{\mathrm{LB}}=-\beta {\left[\rho \left(\mathbf{r}\right)\right]}^{1/3}\frac{{\left(\xi X\right)}^2}{1+3\beta \xi X{sinh}^{-1}\left(\xi X\right)}.$$

(4)

In Equation (4), $\beta =0.01$ and is empirical, while $X=[\nabla \rho ]/{\rho }^{4/3}$. The parameter $\xi$ is a factor that arises in transition from the spin-polarized to spin-unpolarized form [36]. Such a gradient correction approach to the $XC$ functional is built more into the theory and, hence, is less artificial than SIC. We refer to this model of non-local correction as DFT-LB94. It improves the asymptotic behavior of the electron when compared to the exact KS potentials computed from correlated densities. This model potential holds the unitary transformation of all occupied and unoccupied orbitals that are eigenstates of the KS equation with one DFT-LB94 potential. In contrast, the description of DFT-SIC unoccupied orbitals is a bit ambiguous due to the ad hoc nature of the potential in Equation (3).

However, accurate excitation and ionization energies for atoms, molecules, and small clusters [37] are produced using asymptotically correct functionals. Thus, it was anticipated that both DFT-SIC and DFT-LB94 would compete to significantly improve the quality of ground, excited, and continuum spectra in the current study.

2.2. LR-TDDFT Dynamical Response

A time-dependent DFT (TDDFT) approach is used to compute the dynamical response of the clusters to the electromagnetic radiation [38]. The system’s behavior is studied in response to a time-dependent weak external perturbation. The external perturbation z, which represents dipole interaction with linearly polarized light, induces a frequency-dependent complex change in the electron density, $\delta \rho (\mathbf{r};\omega )$. Thus, the linear response of the system can be determined using the density–density response function $\chi$ by

$$\delta \rho (\mathbf{r};\omega )=\int \chi (\mathbf{r},{\mathbf{r}}^{\prime };\omega )z^{\prime }d{\mathbf{r}}^{\prime },$$

(5)

where the full susceptibility $\chi$ includes the electrons’ dynamical correlations. In the auxiliary KS system the same induced density can be equivalently calculated using

$$\delta \rho (\mathbf{r};\omega )=\int {\chi }_0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )\delta V_{eff}({\mathbf{r}}^{\prime };\omega )d{\mathbf{r}}^{\prime },$$

(6)

where $\delta V_{eff}$ includes the external field, plus the induced Hartree and induced $XC$ potentials as follows

$$\delta V_{eff}({\mathbf{r}}^{\prime };\omega )=z^{\prime }+V_{ind}({\mathbf{r}}^{\prime };\omega )$$

(7)

with

$$V_{ind}({\mathbf{r}}^{\prime };\omega )=\int \frac{\delta \rho (\mathbf{r};\omega )}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}+{\left[\frac{\partial V_{XC}}{\partial \rho }\right]}_{\rho ={\rho }_0}\delta \rho ({\mathbf{r}}^{\prime };\omega ).$$

(8)

The response of non-interacting electrons, that is the independent particle (IP) susceptibility, is described by the KS response function ${\chi }_0$, which can be expressed in terms of the ground state KS eigenvalues ${ϵ}_i$ and eigenfunctions ${\varphi }_i$ as

$${\chi }_0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\sum _i{\varphi }_i^{*}\left(\mathbf{r}\right){\varphi }_i\left({\mathbf{r}}^{\prime }\right)G(\mathbf{r},{\mathbf{r}}^{\prime };{ϵ}_i+\omega )+\sum _i{\varphi }_i\left(\mathbf{r}\right){\varphi }_i^{*}\left({\mathbf{r}}^{\prime }\right)G^{*}(\mathbf{r},{\mathbf{r}}^{\prime };{ϵ}_i-\omega ).$$

(9)

Here the index i runs over the occupied states only. For a spherically symmetric atomic cluster, the Green’s function for a parameter E can be expanded in the spherical basis as $G(\mathbf{r},{\mathbf{r}}^{\prime };E)={\sum }_{lm}{G}_{lm}(r,r^{\prime };E)Y_{lm}^{*}(\Omega )Y_{lm}\left({\Omega }^{\prime }\right)$, where the radial component $G_{lm}(r,r^{\prime };E)$ satisfies the radial equation

$$\left(\frac{1}{r^2}\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}-\frac{l(l+1)}{r^2}-V_{DFT}+E\right)G_{lm}(r,r^{\prime };E)=\frac{\delta (r-r^{\prime })}{r^2}$$

(10)

with $G_{lm}(r,r^{\prime };E)=\frac{j_l(r_{<};E)h_l(r_{>};E)}{W{[j_l,h_l]}_{r=c}}$, where $j_l$ and $h_l$ are homogenous solutions of Equation (10) satisfying boundary conditions at $r=0$ and $r=\infty$, respectively, and W is the Wronskian, which is independent of an arbitrary constant c.

A Dyson-like equation for the interacting response function, $\chi$, can be easily derived from Equations (5) and (6) as

$$\chi (\mathbf{r},{\mathbf{r}}^{\prime };\omega )={\chi }_0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )+\int d{\mathbf{r}}^{″}d{\mathbf{r}}^{‴}{\chi }_0(\mathbf{r},{\mathbf{r}}^{″};\omega )\left(\frac{1}{|{\mathbf{r}}^{″}-{\mathbf{r}}^{‴}|}+f_{XC}({\mathbf{r}}^{″},{\mathbf{r}}^{‴};\omega )\right)\chi ({\mathbf{r}}^{‴},{\mathbf{r}}^{\prime };\omega ),$$

(11)

where $f_{XC}$ is the so-called time-dependent $XC$ kernel evaluated with approximate $XC$ functional at ground state density ${\rho }_0$ and expressed as,

$$f_{XC}\left[\rho \right](\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\frac{\delta V_{XC}\left[\rho \right](\mathbf{r};\omega )}{\delta \rho ({\mathbf{r}}^{\prime };\omega )}{|}_{\rho ={\rho }_0}.$$

(12)

Equation (11) can be solved self-consistently for $\chi$ at any desired photon energy. The $\chi$ in Equation (11) can be computed using matrix notation as

$$\chi ={\left(1-\frac{\partial V}{\partial \rho }{\chi }_0\right)}^{-1}{\chi }_0.$$

(13)

where V refers to the ground state potential. Equation (13) can then be solved for $\chi$ using the matrix inversion method [39]. $\delta \rho$ and, hence, $\delta V_{eff}$ can be directly obtained via Equations (5) and (7), respectively.

In LR-TDDFT formalism, the photoionization (PI) cross-section associated with a bound-to-continuum dipole transition is then computed as the sum of independent subshell cross-sections, ${\sigma }_{nl\to k{l}^{\prime }}$, and is given by:

$${\sigma }_{PI}\left(\omega \right)=\sum _{nl}{\sigma }_{nl\to k{l}^{\prime }}\sim \sum _{nl}2(2l+1)|\langle {\psi }_{k{l}^{\prime }}|\delta V_{eff}|{\psi }_{nl}{\rangle |}^2.$$

(14)

The radial component $R_{k{l}^{\prime }}$ of the final continuum wavefunction ${\psi }_{k{l}^{\prime }}$ has the appropriate asymptotic behavior:

$$\underset{r\to \infty }{lim}{R}_{k{l}^{\prime }}\sim \underset{r\to \infty }{lim}[cos\left({\delta }_{l^{\prime }}\right)f_{l^{\prime }}\left(kr\right)+sin\left({\delta }_{l^{\prime }}\right)g_{l^{\prime }}\left(kr\right))] = sin(kr-\frac{l^{\prime }\pi }{2}+\frac{z}{k}ln\left(2kr\right)+{\zeta }_{l^{\prime }}+{\delta }_{l^{\prime }}),$$

(15)

where $f_l$ and $g_l$ represent the regular and irregular spherical Coulomb functions, respectively, and ${\delta }_l$ and ${\zeta }_l=arg\Gamma (l+1-iz/k)$ are, respectively, the short-range and Coulomb phase shifts seen by the ejected electron.

In addition to the external perturbation z, Equation (14) also includes the complex induced field $V_{ind}$ produced by the many-electron interactions. Evidently, the IP level DFT cross-section that disregards correlations is obtained by setting $\delta V_{eff}=z$. This approach makes it simple to compare DFT and TDDFT in order to investigate the role of the many-electron effects during the photoionization process. In this work, we employ the two $XC$ kernels to calculate the PI cross-sections: one by a global averaging procedure, $f_{XC}^{SIC}=\frac{N-1}{N}{V}_{XC}^{LDA}$ with $V_{XC}$ given in Equation (2). We refer to the PI cross-section calculated in this regime as LR-TDDFT-SIC. The other $XC$ kernel with $V_{XC}$ in Equation (2), augmented by Equation (4), yields $f_{XC}^{LB94}$, which in turn is used to evaluate PI cross-sections in LR-TDDFT-LB94.

3. Results and Discussion

3.1. Ground State Structure

In Figure 1, we show the ground state radial potentials of Na${}_{20}$, Na${}_{40}$, and Na${}_{92}$ calculated using DFT-SIC and DFT-LB94. The DFT-SIC curves are obtained by taking an occupancy-weighted average over all the subshells. As the cluster size grows, the potential depth remains roughly unchanged, predominantly since the average density $\rho \left(r\right)$ remains nearly the same. As a consequence, the electron energy levels should become denser with increasing size [40]. This is seen in both SIC and LB94. The shapes of the radial orbitals of the two outer levels (HOMO and HOMO-1), shown in the inset, are almost the same in LB94 and SIC. For a metal cluster, the potential is expected to remain flat in the interior region of delocalized quasi-free electrons while exhibiting strong screening at the edge around R${}_c$ (cluster radius). These characteristics are exhibited by DFT-LB94, while DFT-SIC shows some spatial variations in the interior region. In SIC, these variations are due to the electrons in each orbital experiencing a distinct potential. Additionally, in DFT-SIC, there appear unphysical cusp-like structures below R${}_c$ caused from the radial nodes of occupied subshells [41]. These structures are absent in DFT-LB94, as it is more of an ab initio approach compared to the artificial vacancy-elimination technique in SIC.

The result presented in Figure 1 clearly indicates that in the region outside the radius, the DFT-LB94 potential is somewhat deeper (attractive) and better approximates the $-1/r$ behavior toward the asymptotic region compared to DFT-SIC. In other words, the DFT-LB94 potential has a more accurate asymptotic representation resulting in a slower decay beyond R${}_c$. Such a deeper potential shape allows DFT-LB94 to generate an extended number of virtual unoccupied KS orbitals compared to SIC. On the other hand, it may be noted that the occupied levels in SIC in Na${}_{20}$ are slightly more bound compared to LB94. However, as the cluster size increases, this difference tends to diminish as an indication of a gradual irrelevance of a particular choice of $XC$ with increasing size. Recall that appropriate pseudo-potentials are added to match the ionization thresholds with the experimental values, which is evident from the identical binding energies of HOMO levels obtained in SIC and LB94 in each cluster.

3.2. Total Photoionization Cross-Section

Figure 2 presents the LR-TDDFT photoionization cross-sections for three clusters, Na${}_{20}$, Na${}_{40}$, and Na${}_{92}$, along with the IP results obtained in LR-DFT. TDDFT and DFT calculations are performed using the two $XC$ schemes. The graphs are displayed on a logarithmic scale to emphasize the characteristics at the low-energy region near the ionization threshold. For energies exceeding 20 eV, there are agreements between the TDDFT and DFT results for all three clusters. Furthermore, collective effects disappear leading to discernible oscillations caused by the photoelectron with momentum k being diffracted from the cluster edges [41]. Since this characteristic is related to the geometry of the cluster, which is an IP attribute, they are also observed in the DFT results. Equivalently, the vanishing electron correlations at higher energies is evident from the convergence of oscillations in TDDFT and DFT. The specific choice of an $XC$ functional that is insensitive to diffraction also has no significance at these energies. As the cluster size increases, the frequency of oscillations is also found to increase—a fact that can be explained using a Woods–Saxon model potential. Using this model, the dipole partial photoionization cross-section can be obtained as in [42]

$${\sigma }_{nl\to k{l}^{\prime }}\left(k\right)\propto \frac{e^{-ak}}{{\omega }^{5/2}}(1+cos(2k{R}_c-l^{\prime }\pi )).$$

The decay that follows an exponential pattern is a consequence of the steepness (represented by the parameter “a”) around the edge (refer to Figure 1). Note that potentials being “soft” at the cluster edge is intrinsic in DFT calculations. Based on the equation above, therefore, the state-selected cross-sections ${\sigma }_{nl\to k{l}^{\prime }}$ will exhibit oscillations in k (or, equivalently, in photon energy) with a frequency of 2R${}_c$. Hence, this suggests an increase in oscillations when the cluster size (R${}_c$) increases as seen in the total cross-sections in Figure 2.

When photon energies are in proximity to the ionization threshold and below 10 eV, notable differences emerge between the TDDFT and DFT cross-section profiles as a function of photon energy. These disparities encompass significant enhancements right above the threshold, a host of narrow autoionization resonances, and the occurrence of a correlation minimum around 6–7 eV, all in LR-TDDFT for the three clusters. The tiny and discrete jumps seen in DFT total cross-sections in Figure 2 are due to the opening of new subshell ionization channels.

As noted, the LR-TDDFT cross-sections, regardless of the $XC$ potential chosen, exhibit a strong enhancement in the vicinity of the threshold region for the three clusters compared to LR-DFT. This enhancement can be attributed to the spillover effect of plasmon resonances, which arise from correlated collective electronic motion. Our previous work already characterized this enhancement as the sum of an extension of the surface plasmon resonance from the absorption spectrum and a bulk portion of the volume plasmon resonance [32]. Overall, both the LB94 and SIC approaches exhibit qualitatively similar enhancements. This spillover from discrete to continuum spectra is accompanied by several narrow spikes known as the autoionization resonances, which arise due to the degeneracy between the ionization channels and inner-level single-electron excitation channels.

The phenomenon of plasmon resonance is the result of an in-phase coherence mechanism arising from the coupling of various degenerate ionization channels that are present. In Fano’s framework based on first-order perturbation theory [43], the correlation-modified (LR-TDDFT) matrix element, $M_{nl\to ϵ\lambda }\left(E\right)$, of a dipole ionization channel $nl\to ϵ\lambda$ can be written as

$${\displaystyle M_{n\ell \to ϵ\lambda }\left(E\right)=D_{n\ell \to ϵ\lambda }\left(E\right)+\sum _{n^{\prime }{\ell }^{\prime }\ne n\ell }\int d{E}^{\prime }\frac{\langle {\psi }_{n^{\prime }{\ell }^{\prime }\to {ϵ}^{\prime }{\lambda }^{\prime }}\left(E^{\prime }\right)|\frac{1}{|{\mathbf{r}}_{n\ell }-{\mathbf{r}}_{n^{\prime }{\ell }^{\prime }}|}|{\psi }_{n\ell \to ϵ\lambda }\left(E\right)\rangle }{E-E^{\prime }}\times D_{n^{\prime }{\ell }^{\prime }\to {ϵ}^{\prime }{\lambda }^{\prime }}\left(E^{\prime }\right),}$$

(16)

where $D_{n\ell \to ϵ\lambda }$ refers to the unperturbed (DFT) matrix element. The wavefunctions of the interacting continuum channels are represented by $|\psi \rangle$, and the summation is carried over to all degenerate continuum channels with the exception of the $n\ell$ channel. The significance of electron correlations in enhancing the plasmon resonance spillover in the $n\ell \to ϵ\lambda$ channel can be seen in the second term of Equation (16), which is referred as the interchannel coupling matrix element. ${\mathbf{r}}_{n\ell }$ and ${\mathbf{r}}_{n^{\prime }{\ell }^{\prime }}$ are the spatial co-ordinates of photoelectrons in the interacting continuum channels from initially occupied sub-shells $n\ell$ and $n^{\prime }{\ell }^{\prime }$, respectively. The summation over all subshells exhibits a coherent mixing primarily due to the bound wavefunctions occupying similar spatial regions. This results in a significant increase in the LR-TDDFT cross-section, as shown in Figure 2. However, as illustrated in Figure 2, the DFT predictions that disregard electron correlations do not exhibit such enhancement.

The interpretation of narrow resonances within the Fano framework requires coupling between the bound (excited) and continuum channels, resulting in a modification of the coupling matrix element in Equation (16) as

$$\langle {\psi }_{n^{\prime }{\ell }^{\prime }\to {\eta }^{\prime }{\lambda }^{\prime }}\left(E^{\prime }\right)|\frac{1}{|{\mathbf{r}}_{n\ell }-{\mathbf{r}}_{n^{\prime }{\ell }^{\prime }}|}|{\psi }_{n\ell \to ϵ\lambda }\left(E\right)\rangle$$

where $n^{\prime }{\ell }^{\prime }\to {\eta }^{\prime }{\lambda }^{\prime }$ denote discrete excitation channels. The detailed profiles of these resonant structures are seen in Figure 2 to be sensitive to the scheme of $XC$ functional employed. Resonance positions and shapes differ between SIC and LB94. This happens as a result of the levels of Na${}_{20}$ and Na${}_{40}$ in SIC being energetically slightly deeper compared to their LB94 counterparts, as seen in Figure 1. Additionally, the descriptions of the unoccupied excited states, which depend on the potential’s long-range behavior, differ as well. Since the LB94 potential shows a slightly better long-range behavior, the resonances in this scheme are expected to be somewhat more accurate. The density of the unoccupied levels of both LB94 and SIC being increased as the cluster grows larger (noted in Figure 2) results in a higher number of resonances with a growing size.

3.3. Comparison with Experiments

It is useful, in particular, to assess the effect of the $XC$ schemes on PI below 8 eV, where the collective effect dominates. However, as seen, this region has a complicated spectra because of the presence of autoionization resonances mentioned above. These narrow resonances are usually not present in experimental spectra due to the finite temperature effect experienced by the metal clusters under experimental conditions. In fact, the measured spectrum displays an incoherent mixture of spectra from various satellite configurations driven by the temperature, acquiring, effectively, a specific width [44,45]. This width can camouflage the narrow spikes by smearing them. To simulate this thermal effect, our theoretical data are convoluted with a Gaussian of width 0.4 eV. These are shown for the two clusters in Figure 3, and compared with available measurements. Note, however, that for Na${}_{40}$, experimental data are not available. It can clearly be seen from Figure 3a,b that the smoothed LR-TDDFT cross-sections calculated both in LB94 and SIC show a reasonable agreement with the experimental observations that describe the spillover contribution of the plasmon resonances. Interestingly, SIC produces a slightly better overall agreement, even though its ground state long-range behavior is found to be somewhat less accurate. This suggests a lesser sensitivity of asymptotic properties to the collective behavior.

As can be seen from Figure 3a,b, the peaks of the smoothed LR-TDDFT curve, which must be predominantly induced by the collective volume resonance, are slowly red-shifting as the cluster size increases irrespective of the $XC$ choice. However, the peak position is consistently more red-shifted in LB94 in comparison to SIC for each cluster.

It should be further noticed that the correlation minimum determined in SIC occurs at somewhat higher energies as compared to that in LB94 for all the three systems. As the cluster size increases, this energy separation between the minima from SIC versus LB94 diminishes, while their shapes become narrower. The discontinuities seen around 6 eV in the cross-section results are numerical artifacts, since the convolution is applied until around 6 eV, and therefore, are innocuous. Extending the region of convolution will modify the shape of the correlation minimum that we preferred to avoid.

The present work primarily focuses on the cross-section calculations. It may be interesting to see the effect of XC potentials on the angular distribution parameter since this quantity involves the phase of transition matrix elements. The present work does not make an attempt to study angular distribution because there are no related experimental data on neutral Na clusters. However the angular distribution and the photoelectron spectra measurements on anionic Na${}_n^{-}$ clusters were reported in the past [46,47], and subsequent theoretical studies using RPAE with the jellium model have shown a good agreement with the experiment [48,49].

3.4. Self-Consistent Induced Potential

Important electron correlations are embodied in the dipole matrix element in Equation (14) through the complex induced field, V${}_{ind}$, computed using Equation (8). These correlations play a crucial role for the plasmonic enhancement in the cross-section. The behavior of V${}_{ind}$ delineates the detailed dynamics by visualizing a plasmon as a driven collective-electron oscillator with damping [15,50]. In this model, the real part of V${}_{ind}$, Re(V${}_{ind}$), represents the effective driving field, whereas the imaginary part, Im(V${}_{ind}$), denotes the collective response of the system. Generally, at energies below the resonance peak, Re(V${}_{ind}$) < 0 screens the external field so the electrons can build the collective motion and, simultaneously, the deepening of Im(V${}_{ind}$) favors resonant binding. At the resonance peak, the Re(V${}_{ind}$) becomes zero, making the field irrelevant, where Im(V${}_{ind}$) offers the maximum binding so the collective response perfectly dominates. Above the peak energy, the plasmon decays with increasing Im(V${}_{ind}$) since the field switches to the damping mode (anti-screening) as Re(V${}_{ind}$) > 0. We expect this behavior in our current results but they are dominantly governed by the volume plasmon as the surface plasmon rules only below the ionization threshold.

Figure 4 shows the 3D plots of the real and imaginary components of the radial part of V${}_{ind}$ for all three clusters. Over the plasmon spillover region just above the threshold, Im(V${}_{ind}$) shows a broad well-type shape from the collective dynamics, resulting in a transient attractive force that an emerging photoelectron will feel. In addition, the Re(V${}_{ind}$) switches sign over this range for all three clusters as expected, although the negative Re(V${}_{ind}$) range is shorter, as expected. The details of these results are seen to be sensitive to the form of the $XC$ functional used. Within the graph scale, the magnitudes of the V${}_{ind}$ are slightly higher in LB94 than SIC. This may explain a slightly higher value of the LB94 cross-sections in Figure 3. Moreover, it is observed that the Im(V${}_{ind}$) for SIC are a bit broader. This subsequently sustains the dynamics for a wider energy range to push the correlation minimum to higher energies in SIC, as also seen in Figure 3.

4. Summary and Conclusions

The ground state structures and the photoionization properties of three sodium cluster systems, Na${}_x$ (x = 20, 40, and 92), were studied in a jellium-based LR-TDDFT methodology. Two previously successful schemes, SIC and LB94, of the $XC$ functional implementation were employed. The energy levels obtained from the two schemes differ to some degree, while this difference is found to be at the maximum for Na${}_{20}$. However, the difference reduces as the cluster size increases. The occupied level distributions remain roughly the same in the two schemes for each cluster. Both of the $XC$ schemes show features such as the above-threshold enhancement from the plasmon resonant spillover, the narrow single-electron autoionizing resonances, the appearance of the correlation minimum, and the diffractive oscillations at higher energies. However, the $XC$ functionals compete in producing the detailed characteristics of these features. The plasmon-enhanced spillover spectra in both $XC$ formalisms show reasonable agreement with the available experiments. The results further indicate that the cluster size alters the ground state and ionization properties fairly monotonically, but the incongruence between SIC and LB94 tends to close as the size grows.

Even though the quantitative details are found to be sensitive to the particular $XC$ scheme chosen, the significant qualitative similarities provide confidence on the accuracy of a DFT-level description of a metal cluster’s static and dynamical properties. We hope the present work will motivate new experiments in the cluster science, particularly using the photoelectron spectroscopic techniques, which can provide further information in optimizing computational methods. This will help in extending calculations to address, for instance, the noble metal clusters and even to technologically relevant complexes with metal clusters as seed components. The current study may also motivate molecular-level calculations, going beyond the jellium model, to test the outcomes specifically and to increase the accuracy in general, although describing the continuum in such a frame will be a steep challenge to overcome.