Study of high-order harmonic generation in xenon based on time-dependent density-functional theory

High-order harmonic generation (HHG) is one of the most well-known phenomena of strong-field physics, which is of interest for numerous applications [1, 2]. HHG serves as a unique tool for producing subfemtosecond and attosecond highly coherent light pulses in the extreme ultraviolet (XUV) range. The high intensity of generated XUV pulses and their precise synchronization with the laser field allows pump-probe studies of ultrafast phenomena in quantum systems from atoms and molecules to periodic solid structures (see [3, 4] and references therein). Moreover, XUV pulses obtained within HHG can be used for far-field imaging with nanometer resolution [5], frequency-stabilized seeding for free-electron lasers [6] and controlling magnetic dynamics [7].

The physical mechanism of HHG is qualitatively explained in the framework of a three-step scenario [8]. At the first step, atom or molecule is ionized by an intense field, at the second step, liberated electron is accelerated, and at the third step, electron returns back and recombines with the emission of high-energy photon in a wide energy range. While at the second step, the decisive factor is the interaction of moving electron with the laser field, the type and state of the target affect the ionization probability and the recombination cross section [9–20]. Thus, the multielectron dynamics associated with the motion of various electrons in an atom or molecule can significantly affect both the form of HHG spectra and the optimal generation regimes [12, 14, 21–29].

A typical example of the large influence of multielectron effects at the recombination stage is the giant enhancement of HHG yield in xenon, associated with the interaction of the photoelectron, originating from the external 5p₀ orbital, with the inner electrons of the atom [14, 23–26, 30, 31]. This enhancement was first predicted in [30] based on the analytical parametrization of the HHG spectra within the so-called electron wave packet (EWP), which describes the ionization and propagation of photoelectron in the continuum, and a differential photorecombination cross section (PRCS), which is related to a differential photoionization cross section (PICS) according to the principle of the detailed balance [32]. PICS from the outer 5p subshell of Xe atom has a broad maximum near photon energy of 100 eV due to the giant resonance in the PICS from the inner 4d subshell and interchannel coupling via electron–electron interaction [33, 34]. Experimental measurements of HHG spectra in Xe using near-IR driving wavelengths confirm the existence of a giant (more than an order of magnitude) harmonic yield enhancement near 100 eV, which qualitatively agrees with the result of analytical parameterization [10, 14, 26, 30, 31]. The giant enhancement in the HHG spectrum of Xe atom was numerically studied in references [24, 26] by applying time-dependent configuration-interaction singles approach and in reference [25] by using time-dependent R-matrix theory. In these studies, the dynamical interaction between orbitals on the outer 5th shell and 4d subshell was taken into account, while the other orbitals were freezed. Although the calculated HHG spectra's coincidence with results of analytical parametrization and experiment is not perfect, obtained HHG spectra contain a broad maximum in the vicinity of 100 eV, which disappears if the interaction between the orbitals is excluded from calculations. Moreover, all five 4d_m orbitals (m = 0, ±1, ±2) on 4d subshell contribute to the HHG, and it is not sufficient to take into account only the 4d₀ orbital aligned in the direction of the laser polarization [24]. However, due to the rapid increase in the required computational resources with the number of active orbitals, the role of other orbitals on the 4th shell in HHG has not been studied previously.

An alternative approach that may be applicable to investigate the enhancement of HHG yield in Xe is the time-dependent density-functional theory (TDDFT) [35, 36]. The heart of the TDDFT are the time-dependent Kohn–Sham (TDKS) equations for atomic orbitals, which include the interaction of electrons with the laser pulse, nucleus, and electron–electron exchange–correlation interaction. This approach describes with high accuracy giant resonance in the PICSs of various Xe and Xe-like ions [37–40]. However, this approach has not been used previously to study HHG in Xe due to the high complexity of calculations. The main difficulties are caused by a high probability of recolliding electron to move over long distances and acquire high energies, the very sharp increase of the potential energy in the vicinity of the ion, the divergences arising at partial freezing of individual shells of an atom.

Recently, we developed the numerical code for TDKS equations solution in a spherical coordinate system that is based on the careful discretization of the numerical grid, careful selection of calculation region size, and parallelization between the nodes of a computing cluster [20]. This code allowed to simulate HHG in argon taking into account all electrons and demonstrate that a significant multielectron effect in HHG at high laser-pulse intensity is the laser-induced polarization of the atom itself [20]. In this work, we use this numerical code to carry out numerical simulation of the HHG in Xe, taking into account the dynamics of all subshells on the 4th and 5th atomic shells. We show that at the enhancement region and at higher frequencies HHG is influenced not only by the 4d subshell, but also the other, 4p and 4s subshells on the 4th shell. The results are explained based on the dynamics of electron–electron potential induced by the recombining 5p₀ wave packet and the interaction between electrons on the 4th shell. We also perform the analysis of the Xe atom interaction with monochromatic XUV field and demonstrate a high linear response of 4p and 4s subshells.

The paper is organized as follows. In section 2, we discuss the TDKS equations and details of the numerical method. Section 3 is devoted to study the dynamics of orbitals of Xe atom in monochromatic XUV field. In section 4 we calculate HHG spectrum in the infrared field and compare it with the results of experiment and the analytical parametrization of the HHG spectrum, which provides high accuracy of EWP calculation at frequencies below the cutoff frequency [41]. Contributions of different orbitals and subshells to the total dipole acceleration are studied in section 4.1. Our main results are summarized in section 5. In the appendix we present the details of the analytical parametrization of the HHG spectrum. Atomic units (a.u.) are used throughout the paper unless specified otherwise.

2.1. Initial electronic structure of xenon atom

The initial condition for the TDKS equations corresponds to the unperturbed Xe atom in the ground state. It is found within stationary density functional theory, which expresses the state of a multielectron atom in terms of stationary Kohn–Sham (KS) orbitals, ${\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)$. The unperturbed KS orbitals satisfy the system of stationary KS equations:

$$\begin{equation}{\hat{H}}_{0}{\psi }_{j}^{\left(0\right)}={E}_{j}{\psi }_{j}^{\left(0\right)},\quad j=1,2,\dots ,N/2,\end{equation} \tag{ 2.1 }$$

$$\begin{equation}{\hat{H}}_{0}=-\frac{1}{2}{\nabla }^{2}-\frac{N}{r}+{V}_{\text{ee}}\left[{\rho }_{0}\left(\mathbf{r}\right)\right],\end{equation} \tag{ 2.2 }$$

where N is the total number of electrons in an atom (for Xe, N = 54), V_ee[ρ₀] = V_H[ρ₀] + V_xc[ρ₀] is the electron–electron interaction potential, which is the functional of the electron density ⁵

$$\begin{equation}{\rho }_{0}\left(\mathbf{r}\right)=2\sum\limits _{j=1}^{N/2}{\left\vert {\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)\right\vert }^{2}.\end{equation} \tag{ 2.3 }$$

The electron–electron interaction potential is presented by a sum of Hartree potential describing the electron repulsion in the mean-field approximation,

$$\begin{equation}{V}_{\mathrm{H}}\left[{\rho }_{0}\right]\left(\mathbf{r}\right)=\int \frac{{\rho }_{0}\left({\mathbf{r}}^{\prime }\right)}{\left\vert \mathbf{r}-{\mathbf{r}}^{\prime }\right\vert }\mathrm{d}{\mathbf{r}}^{\prime },\end{equation} \tag{ 2.4 }$$

and the exchange–correlation potential V_xc[ρ₀], which is considered in the spin-unpolarized LB94 approximation [42]. Since all shells of Xe atom are fully occupied, both ρ₀(r) and V_ee[ρ₀] are spherically symmetric, so that a particular solution of the equation (2.1) can be characterized by the principal quantum number (n), angular momentum (l) and its magnetic projection (m) on the quantization axis z:

$$\begin{equation}{\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)={r}^{-1}{{\Psi}}_{nl}^{\left(0\right)}\left(r\right){Y}_{lm}\left(\theta ,\varphi \right),\end{equation} \tag{ 2.5 }$$

where ${{\Psi}}_{nl}^{\left(0\right)}\left(r\right)$ is the radial part of KS orbital, Y_lm (θ, φ) is the spherical harmonic, angles θ and φ are polar and azimuthal angles of the spherical coordinate system. Principal quantum number $n=\bar{1,5}$ corresponds to the number of electron shell; l corresponds to different subshells and $l=\bar{0,\left(n-1\right)}$ for n ⩽ 3, l = {0, 1, 2} for n = 4, and l = {0, 1} for n = 5; m changes from −l to l for any subshell. The energy levels are degenerate in m, so that the orbital wavefunctions ${\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)$ for fixed subshell (n, l) correspond to the same binding energy. For the definiteness we assume that KS orbital index j increases with n, l and m, i.e.

$$\begin{equation}j=\begin{cases}n\left(n-1\right)\left(2n-1\right)/6+m+l+{l}^{2}+1\quad \hfill & \quad n{\leqslant}4,\enspace l{< }\mathrm{min}\left(n,3\right),\hfill \\ m+l+{l}^{2}+24\quad \hfill & \quad n=5,\enspace l{< }2.\hfill \end{cases}\end{equation} \tag{ 2.6 }$$

The corresponding ground state configuration, taking into account the presence of two electrons with different spins in each state (n, l, m), is written as [Ar]3d¹⁰4s²4p⁶4d¹⁰5s²5p⁶, where [Ar] = 1s²2s²2p⁶3s²3p⁶ is the electronic configuration of Ar atom.

The radial parts of KS orbitals Ψ_nl⁽⁰⁾ (r) are found within the imaginary time propagation method [43] (see the numerical scheme of time propagation in section 2.2). In order to resolve oscillations of wave functions with the scale ∼1/N near the nucleus, we use a nonuniform spatial grid, which has a higher density of nodes near the nucleus. The radial nodes of the spatial grid are specified as

$$\begin{equation*}{r}_{k}=k{\Delta}r+\left(\delta r/{\Delta}r-1\right){r}_{\alpha }\enspace \mathrm{tanh}\left(k{\Delta}r/{r}_{\alpha }\right),\end{equation*}$$

where k is an integer, δr = 10⁻⁴ a.u. is tiny step of the radial grid, which is realized near the nucleus; Δr = 0.1 a.u. is the radial step for large distances; r_α = 20 a.u. is the scale for changing spatial step from δr to Δr. The size of radial grid was limited by 200 a.u.

The binding energies of KS orbitals are derived from (2.1). They are listed in table 1 along with experimental binding energies, which resolve spin–orbit splitting of p and d energy levels [44, 45]. The deviation of the numerical and experimental binding energies do not exceed 11%. In particular, for the outer 5p subshell, the accuracy of agreement with the experiment is 5%, for 5s, 4p, and 4d subshells, it is 1%, and for the 4s orbital, it is 11%. Thus, the used LB94 exchange–correlation approximation is generally well suited for describing the electronic configuration of the Xe atom.

Table 1. Absolute values of orbital binding energies of Xe. (A) Our DFT calculation with LB94 exchange–correlation potential [42]. (B) Experiment [44, 45].

Orbital	Energy (eV)
	(A) DFT	(B) Experiment
1s	33 068.8	34 561
2s	5032.4	5453
2p	4741.3	5107 (3p1/2), 4786 (3p3/2)
3s	1039.8	1148.7
3p	915.2	1002.1 (3p1/2), 940.6 (3p3/2)
3d	682.6	689 (3d3/2), 676.4 (3d5/2)
4s	189.4	213.2
4p	145.5	146.7 (4p1/2), 145.5 (4p3/2)
4d	69.9	69.5 (4d3/2), 67.5 (4d5/2)
5s	23	23.3
5p	12.8	13.4 (5p1/2), 12.1 (5p3/2)

2.2. Time-dependent Kohn–Sham equations

We consider the interaction of a linearly polarized laser pulse with an atomic target in the dipole approximation. For this case, TDKS equations are

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}{\psi }_{j}\left(\mathbf{r},t\right)=\hat{H}{\psi }_{j}\left(\mathbf{r},t\right),\end{equation} \tag{ 2.7 }$$

$$\begin{equation}\hat{H}=-\frac{{\nabla }^{2}}{2}-\frac{N}{r}+\mathbf{r}\cdot \mathbf{F}\left(t\right)+{V}_{\text{ee}}\left[\rho \left(\mathbf{r},t\right)\right],\end{equation} \tag{ 2.8 }$$

where $\mathbf{F}\left(t\right)=\hat{\mathbf{z}}F\left(t\right)$ is the electric field of a short laser pulse, ψ_j (r, t) is the wave function of the TDKS orbital determining the time-dependent electron density

$$\begin{equation}\rho \left(\mathbf{r},t\right)=2\sum\limits _{j=1}^{N/2}{\left\vert {\psi }_{j}\left(\mathbf{r},t\right)\right\vert }^{2}.\end{equation} \tag{ 2.9 }$$

The TDKS equations are solved with initial condition

$$\begin{equation}{\psi }_{j}\left(\mathbf{r},t=0\right)={\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right),\end{equation} \tag{ 2.10 }$$

which we have discussed in the previous subsection.

Interaction of an intense laser pulse with an atomic target induces a time-dependent dipole moment D(t). For linearly polarized field, the dipole moment is directed along z axis, $\mathbf{D}\left(t\right)=\hat{\mathbf{z}}D\left(t\right)$ and is expressed in terms of density distribution:

$$\begin{equation}D\left(t\right)=\int z\rho \left(\mathbf{r},t\right)\mathrm{d}\mathbf{r}=\sum\limits _{j}{D}_{j}\left(t\right),\end{equation} \tag{ 2.11 }$$

$$\begin{equation}{D}_{j}\left(t\right)=2\int z\vert {\psi }_{j}\left(\mathbf{r},t\right){\vert }^{2}\enspace \mathrm{d}\mathbf{r}.\end{equation} \tag{ 2.12 }$$

The properties of the radiation generated by atom subjected into intense laser field are determined by the dipole acceleration a(t):

$$\begin{equation}a\left(t\right)={\mathrm{d}}^{2}D\left(t\right)/\mathrm{d}{t}^{2}=\sum\limits _{j}{a}_{j}\left(t\right),\end{equation} \tag{ 2.13 }$$

$$\begin{equation}{a}_{j}\left(t\right)={\mathrm{d}}^{2}{D}_{j}\left(t\right)/\mathrm{d}{t}^{2}=-2F\left(t\right)-2\int \vert {\psi }_{j}{\vert }^{2}\frac{\partial }{\partial z}\left({V}_{\text{ee}}-\frac{N}{r}\right)\mathrm{d}\mathbf{r},\end{equation} \tag{ 2.14 }$$

where a_j (t) is the partial dipole acceleration of the KS orbital having index j. The spectral density of the dipole acceleration is given by squared Fourier transform of a(t),

$$\begin{equation}R\left[a\right]\left(\omega \right)={\left\vert {\int }_{-\infty }^{\infty }a\left(t\right){\mathrm{e}}^{\mathrm{i}\omega t}\enspace \mathrm{d}t\right\vert }^{2}.\end{equation} \tag{ 2.15 }$$

The general formulation of the density functional method allows one to study the contributions of different KS orbitals, subshells, and whole shells of an atom to the total dipole acceleration; thereby, TDKS equations provide access to the inner details of the atomic response in the external field. Corresponding contributions from shell or subshell, respectively, are given by proper summation in j:

$$\begin{equation}{\hat{a}}_{n}\left(t\right)=\sum\limits _{\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{l}\mathrm{l}\enspace n}{a}_{j}\left(t\right),\hspace{25.0pt}{\bar{a}}_{nl}\left(t\right)=\sum\limits _{\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{l}\mathrm{l}\enspace \left(n,l\right)}{a}_{j}\left(t\right),\end{equation} \tag{ 2.16 }$$

where the range of j determines the contribution for whole shell n or either subshell (n, l). The spectral density of partial contributions of KS orbitals, subshells, and shells are calculated using equation (2.15) as R[a_j ](ω), $R\left[{\bar{a}}_{nl}\right]\left(\omega \right)$, and $R\left[{\hat{a}}_{n}\right]\left(\omega \right)$, respectively.

The algorithm for the numerical solution of TDKS equation (2.7) is the same as that used by us in reference [20]. Taking into account that for linearly polarized laser pulse the magnetic quantum number is conserved for each KS orbital, the KS orbitals are expanded in spherical harmonics [43, 46]:

$$\begin{equation}{\psi }_{j}\left(\mathbf{r},t\right)={r}^{-1}\sum\limits _{{l}^{\prime }=0}^{{l}_{\mathrm{max}}}{{\Psi}}_{j{l}^{\prime }}\left(r,t\right){Y}_{{l}^{\prime }m}\left(\theta ,\varphi \right).\end{equation} \tag{ 2.17 }$$

Each term corresponds to the orbital momentum l', and Ψ_jl'(r, t) is the radial part satisfying the initial condition ${{\Psi}}_{j{l}^{\prime }}\left(r,0\right)={\delta }_{{l}^{\prime }l}{{\Psi}}_{nl}^{\left(0\right)}\left(r\right)$. Here, δ_l'l is the Kronecker symbol, and n, l, m are quantum numbers that define KS orbital number j by (2.6). Since |ψ_j (r, t)|² does not depend on the azimuthal angle, the electron–electron interaction potential can be expanded in the Legendre polynomial, P_k (cos θ), series. For considered laser pulse parameters, this expansion may be limited by the first three terms:

$$\begin{equation}{V}_{\text{ee}}\left[\rho \left(\mathbf{r},t\right)\right]\approx \sum\limits _{k=0}^{2}{V}_{\text{ee}}^{\left(k\right)}\left(r,t\right){P}_{k}\left(\mathrm{cos}\theta \right).\end{equation} \tag{ 2.18 }$$

Using expansion (2.17) and approximation (2.18) we reduce equation (2.7) to the system of coupled equations:

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}{\mathbf{\Psi }}_{j}\left(r,t\right)=\left({\hat{\mathbf{R}}}_{j}+{\hat{\mathbf{W}}}_{j}\right){\mathbf{\Psi }}_{j}\left(r,t\right),\end{equation} \tag{ 2.19 }$$

where Ψ_j is the column-vector whose l'th element is Ψ_jl'(r, t). Elements of the diagonal matrix ${\hat{\mathbf{R}}}_{j}$ are

$$\begin{equation}{\left[{\hat{\mathbf{R}}}_{j}\right]}_{{l}^{\prime }{l}^{\prime }}=-\frac{1}{2}\frac{{\partial }^{2}}{\partial {r}^{2}}+\frac{{l}^{\prime }\left({l}^{\prime }+1\right)}{2{r}^{2}}-\frac{N}{r}+{V}_{\text{ee}}^{\left(0\right)}+\frac{{l}^{\prime }\left({l}^{\prime }+1\right)-3{m}^{2}}{\left(2{l}^{\prime }+1\right)\left(2{l}^{\prime }+3\right)}{V}_{\text{ee}}^{\left(2\right)},\end{equation} \tag{ 2.20 }$$

and elements of the pentadiagonal matrix ${\hat{\mathbf{W}}}_{j}$ are determined by

$$\begin{equation}{\left[{\hat{\mathbf{W}}}_{j}\right]}_{k{l}^{\prime }}=\begin{cases}_{{l}^{\prime }-2,m}{V}_{\text{ee}}^{\left(2\right)},\quad \hfill & \quad k={l}^{\prime }-2\hfill \\ {c}_{{l}^{\prime }-1,m}\left[rF\left(t\right)+{V}_{\text{ee}}^{\left(1\right)}\right],\quad \hfill & \quad k={l}^{\prime }-1\hfill \\ {c}_{{l}^{\prime }m}\left[rF\left(t\right)+{V}_{\text{ee}}^{\left(1\right)}\right],\quad \hfill & \quad k={l}^{\prime }+1\hfill \\ {d}_{{l}^{\prime }m}{V}_{\text{ee}}^{\left(2\right)},\quad \hfill & \quad k={l}^{\prime }+2\hfill \\ 0\quad \hfill & \quad \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{cases}\end{equation} \tag{ 2.21 }$$

where

$$\begin{equation}{c}_{lm}=\sqrt{\frac{{\left(l+1\right)}^{2}-{m}^{2}}{\left(2l+1\right)\left(2l+3\right)}},\hspace{25.0pt}{d}_{lm}=\frac{3}{2}{c}_{lm}\sqrt{\frac{{\left(l+2\right)}^{2}-{m}^{2}}{\left(2l+3\right)\left(2l+5\right)}}.\end{equation} \tag{ 2.22 }$$

The propagation of Ψ_j (r, t) over time step Δt is performed using three-split-operator symmetric decomposition [47]:

$$\begin{equation}{\mathbf{\Psi }}_{j}\left(r,t+{\Delta}t\right)\approx {\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{W}}}_{j}\left(t\right)/2}\enspace {\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{R}}}_{j}\left(t\right)}\enspace {\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{W}}}_{j}\left(t\right)/2}{\mathbf{\Psi }}_{j}\left(r,t\right).\end{equation} \tag{ 2.23 }$$

First and third exponential operators are applied through the diagonalization of matrix ${\hat{\mathbf{W}}}_{j}$ within unitary transformation ${\hat{\mathbf{W}}}_{j}={\hat{\mathbf{S}}}_{j}{\hat{\mathbf{W}}}_{j,\mathrm{diag}}{\hat{\mathbf{S}}}_{j}^{{\dagger}}$, where ${\hat{\mathbf{W}}}_{j,\mathrm{diag}}$ is the diagonal matrix and unitary matrix ${\hat{\mathbf{S}}}_{j}$ is known analytically:

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{W}}}_{j}/2}={\hat{\mathbf{S}}}_{j}\enspace {\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{W}}}_{j,\mathrm{diag}}/2}{\hat{\mathbf{S}}}_{j}^{{\dagger}}.\end{equation} \tag{ 2.24 }$$

The second exponential operator is approximated by Crank–Nicolson method:

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}{\Delta}t{\hat{\mathbf{R}}}_{j}}\approx {\left(1+\mathrm{i}{\hat{\mathbf{R}}}_{j}{\Delta}t/2\right)}^{-1}\left(1-\mathrm{i}{\hat{\mathbf{R}}}_{j}{\Delta}t/2\right).\end{equation} \tag{ 2.25 }$$

The second derivative in ${\hat{\mathbf{R}}}_{j}$ is calculated by using the Numerov approximation [46].

The maximum orbital momentum in expansion of wave functions in spherical harmonics is l_max = 64 for calculations in section 3 and l_max = 512 for calculations in section 4. The orbitals on 1–3 shells are frozen, i.e. ${\psi }_{j}\left(\mathbf{r},t\right)\equiv {\mathrm{e}}^{-\mathrm{i}{E}_{j}^{\left(0\right)}t}{\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)$ for j ⩽ 14. Freezing of internal orbitals makes it possible to use a relatively large time step Δt = 0.02 a.u. without losing accuracy in numerical calculations. The absorbing layer is modeled by multi-hump imaginary potentials [48] with total absorption layer width r_abs = 50 a.u., providing a wide range of effectively absorbed de Broglie wavelengths.

In this section, we consider the interaction of Xe atom with a perturbative monochromatic XUV field. Within the numerical solution of TDKS equation (2.7), we calculate the total (angle-integrated) PICS and study the response of individual subshells from the 4th and 5th shells.

The XUV field F(t) is parameterized as F(t) = F₀ f(t)cos(ω₀ t), where ω₀ is the carrier frequency, and f(t) is trapezoidal envelope with duration τ_p = 5 fs for flat top part and τ₀ = 0.5 fs for tuning on/off part. The intensity of the XUV pulse is fixed, I = 2 × 10¹² W cm⁻², while frequency ω₀ is considered in the range 20–220 eV. For such parameters we can neglect the contribution from the two-photon ionization.

The PICS ${\sigma }_{nl}^{\text{ion}}\left({\omega }_{0}\right)$ is found based on the sum of the averaged ionization probabilities per unit time from individual KS orbitals:

$$\begin{equation}{\sigma }_{nl}^{\text{ion}}\left({\omega }_{0}\right)=\frac{{\omega }_{0}}{I}\sum\limits _{\text{subshell}\;\left(n,l\right)}{\bar{w}}_{j}\left({\omega }_{0}\right),\end{equation} \tag{ 3.1 }$$

$$\begin{equation}{\bar{w}}_{j}\left({\omega }_{0}\right)=\frac{2}{{\tau }_{p}}\left(1-\underset{r{< }{r}_{0}}{\int }\vert {\psi }_{j}\left(t={t}_{f}\right){\vert }^{2}\enspace \mathrm{d}\mathbf{r}\right),\end{equation} \tag{ 3.2 }$$

where, t_f = τ_p + 2τ₀ is the ending time of the XUV pulse, r₀ = 10 a.u. determines the radius of the sphere, whose edge separates bounded (r < r₀) and free electrons (r > r₀). The used value of r₀ was chosen based on the stability of ${\sigma }_{nl}^{\text{ion}}$ when changing r₀ near 10 a.u. for all considered n, l in a wide range of ω₀. Figure 1(a) shows PICS ${\sigma }_{nl}^{\text{ion}}\left({\omega }_{0}\right)$ for all subshells on 4th and 5th shells. PICSs for 5p, 5s, and 4d subshells contain wide maxima near ω₀ ≈ 100 eV. The corresponding maximum in 4d PICS is more than an order of magnitude higher than that for 5p and 5s subshells. Note that the calculated PICSs are in excellent agreement with the experimental PICS from reference [49], which confirms the high accuracy of the TDKS equation (2.7) with used electron–electron interaction potential for studying phenomena caused by electron–electron interactions in Xe atom. PICSs for 4s and 4p subshells are negligible in the region of the giant resonance 60 eV < ω < 130 eV due to insufficient photon energy for ionization (the threshold of the photoelectric effect for 4p subshell is |E_4p | = 145.5 eV, for 4s, |E_4s | = 189.4 eV, see table 1).

Zoom In Zoom Out Reset image size

In the central mean-field model, the shape resonance in the 4d PICS originates from a sharp resonance-like increase in the matrix element of the transition from the 4d state to the continuum ɛf state. This results from the double-well structure of the effective potential created by the ion field and centrifugal force. The barrier between these two wells suppresses the overlapping of continuum ɛf and 4d bound state wave functions. This overlapping gradually increases with increasing ɛ and is maximized for those energies, which are located near the maximum of this barrier [50–54]. As a result, the matrix element of the 4d → ɛf transition increases sharply. The further increasing photoelectron energy in f-channel leads to a decreasing transition matrix element due to the nodal structure of 4d orbital, thereby forming the shape of the giant resonance in the 4d PICS.

Our TDKS calculations of PICSs with the frozen potential for electron–electron interaction (V_ee(r, t) ≡ V_ee[ρ₀(r)]) show: (i) the maximum of PICS for 4d subshell is located at ω₀ ≈ 80 eV; (ii) PICSs for 5p and 5s subshells do not show up maxima (see dashed lines in figure 1(a)). The account of the time-dependent dynamics into potential V_ee(r, t) shifts the maximum of the 4d PICS and makes the shape of resonance wider. Moreover, due to the energy exchange with the external shell, PICSs for 5p and 5s subshells also get a strong maximum around 90 eV (see solid lines in figure 1(a)). All these issues are found in agreement with previous results in reference [40].

The broadening and shift of the giant resonance in 4d PICS is primarily the consequence of the 'collectivization' of ten 4d electrons via the interchannel coupling [33, 40, 51, 53, 54]. However, when the Xe atom is exposed to an external field at the frequency ω₀  ≳ 100 eV, a linear response (i.e. induced dipole moment) is observed not only for 4d but also for 4p and 4s subshells. This can be seen from the calculation of spectral densities at the carrier frequency ω₀ of the dipole acceleration for individual subshells, $R\left[{\bar{a}}_{nl}\right]\left({\omega }_{0}\right)$, and for the 4th shell, $R\left[{\hat{a}}_{n}\right]\left({\omega }_{0}\right)$, where n = 4. Freezing of the electron–electron potential leads to a decrease in the 4p and 4s orbitals response in the range ω₀ = 100–130 eV by an order of magnitude or even higher (see figure 1(c)). Thus, the electron–electron interaction with the 4d subshell causes the enhancement of the response of 4p and 4s subshells. This strong interchannel coupling is realized due to similar spatial distribution of all orbitals on the 4th shell (see figure 2). Indeed, according to the perturbation theory for the TDKS equations [37, 55], the coupling matrix element between different KS orbitals is maximized for similar spatial distributions of their wavefunctions. Figure 2 shows that the wide maxima of the radial probability distributions of KS orbitals |Ψ_nl⁽⁰⁾ |²(r) for 4d, 4p and 4s subshells are located at r ∼ 0.68–0.75 a.u., and the characteristic scales of the decreasing from the maximum for 4d, 4p and 4s wave functions are approximately the same. This promotes a laser photon energy transfer from the 4d subshell to orbitals on 4p and 4s subshells through the electron–electron interaction. This physical channel, along with the interaction between electrons on the 4d subshell, can also be one of the reasons for the shift and expansion of the giant resonance in the PICS for the 4d subshell.

Zoom In Zoom Out Reset image size

It should be noted that the calculated responses of individual subshells of the Xe atom in the XUV field contain narrow resonance maxima seen in figure 1(b). The corresponding frequencies equal to the difference in the binding energies of the given subshell and the other subshells having angular momentum differing by one. This resonance maxima originate from the transition between the initial state (n, l, m) of each KS orbital ψ_j to the state (n', l', m) of unperturbed Hamiltonian under the action of external XUV field with a resonant frequency ω₀ = |E_j' − E_j |, where E_j' is the energy of the (n', l', m) level and l' = l ± 1. ⁶ As a result, the dipole moment of jth KS orbital, D_j (t), contains a cross-term that oscillates with this frequency. However, the dipole moment of j'th KS orbital, D_j'(t), also contain oscillations at the same resonant frequency and the same amplitude but the opposite sign, which is due to the back transition (n', l', m) → (n, l, m) in the j'th KS orbital ψ_j'. As a result, all the 'artificial' resonances of KS orbitals cancel each other in the total dipole moment. Indeed, as is seen from figure 1(b), the responses of 4d and 4p subshells have narrow resonance maxima at ω₀ = |E_4d − E_4p | ≈ 75.6 eV, but there is no any enhancement in the total response of the 4th shell at this resonance frequency.

The nonlinear response of the xenon atom appears at high harmonics of the IR frequency when the atom is subjected into an intense IR field. This section considers the HHG in Xe subjected to the two-cycle IR-field having intensity I = 1.9 × 10¹⁴ W cm⁻² and wavelength λ = 1.8 μm. The selected parameters correspond to the experiment [14]. In the numerical calculations, we parameterize the electric component of the laser pulse in terms of vector potential:

$$\begin{equation}F\left(t\right)=-\frac{1}{c}\frac{\partial A}{\partial t},\hspace{25.0pt}A\left(t\right)=-\frac{c{F}_{0}}{{\omega }_{0}}f\left(t\right)\mathrm{sin}\left({\omega }_{0}t\right),\end{equation} \tag{ 4.1 }$$

where ω₀ = 2πc/λ is the carrier frequency, c is the speed of light, F₀ is the strength of a laser pulse, whose magnitude determines the intensity of a laser pulse as $I={F}_{0}^{2}{I}_{\mathrm{a}}$, where F₀ is normalized to the atomic electric-field unit and I_a = 3.51 × 10¹⁶ W cm⁻². The envelope of the vector potential is f(t) = sin²(πt/τ) for 0 < t < τ and f(t) = 0 otherwise. The intensity full-width at half maximum duration τ_p = τ[2 arccos(2^−1/4)/π] = 12 fs corresponds to two field periods at the carrier frequency.

In figure 3(a) we show a comparison of HHG spectra R[a](ω) calculated within TDDFT (blue solid line), analytical parameterization (red solid line) suggested in references [10, 56, 57], and experimental result [14] (black line). ⁷ The analytical parameterization of HHG yield consists in the separation of the atomic and laser parts (EWP):

$$\begin{equation}R\left(\omega \right)=\mathcal{W}\left(\mathcal{E}\right){\sigma }_{\text{rec}}\left(\mathcal{E},0\right),\quad \end{equation} \tag{ 4.2 }$$

where $\mathcal{W}\left(\mathcal{E}\right)$ is the EWP, ${\sigma }_{\text{rec}}\left(\mathcal{E},0\right)$ is the differential PRCS of electron with kinetic energy $\mathcal{E}=\omega -{I}_{p}$, where I_p is the ionization potential (see more details in appendix A). All three HHG spectra are found in good agreement with each other. The spectral intensity (averaged over 2ω₀ range) contains a giant enhancement region with a maximum at ω ≃ 100 eV. In the frequency range below the giant enhancement, the spectrum contains a minimum at ω ≈ 60 eV and a narrow near-threshold enhancement region near ω ≈ 20 eV. For harmonic energies above the giant resonance maximum, the numerical HHG spectrum has a pronounced minimum at ω ≈ 160 eV and gradually increases with ω up to the cutoff at ω ≈ 185 eV. Experimental harmonic yield rapidly drops off at frequencies above the giant enhancement region, which can be caused by the propagation macroscopic effects associated with phase mismatch, defocusing, and others in the experiment [58–61]. Also, our calculations do not take into account the spatial intensity distribution of the laser pulse, which can be also one of the reasons for deviation of HHG spectrum from the experimental results.

Zoom In Zoom Out Reset image size

To ensure that the high harmonic generation in Xe is provided within the standard three-step scenario (and impact ionization does not play an important role), we calculate the Gabor transform of the dipole acceleration found from TDDFT:

$$\begin{equation}S\left(\omega ,t\right)={\int }_{-\infty }^{\infty }a\left(t-{t}^{\prime }\right){\mathrm{e}}^{-\mathrm{i}\omega {t}^{\prime }-\left(2\enspace \mathrm{ln}\enspace 2\right){t}^{\prime 2}/{\tau }_{0}^{2}}\enspace \mathrm{d}{t}^{\prime },\end{equation} \tag{ 4.3 }$$

where τ₀ = 0.2 fs. Figure 3(b) shows the time-frequency distribution |S(ω, t)|² and the classical dependences of the generation frequency for the trajectories with one and higher-order returns [63] (dotted and dashed lines, respectively) on the return time moment. The generation frequency was found as ω_c = E_kin + I_p , where I_p = 12.1 eV is the ionization potential of Xe atom, E_kin is the gained kinetic energy calculated using classical Newton equation (see appendix A). Over the entire frequency range ω > 20 eV, high harmonics are generated along classical trajectories supporting three-step scenario of HHG. An increase in the spectral intensity is observed along all trajectories near ω ≃ 100 eV, which indicates the high amplitude of recombination in the region of giant HHG enhancement. In addition, the spectral intensity is locally increased in the vicinity of caustics corresponding to the merging of 'short' and 'long' closed electrons trajectories (classified in terms of return time) [64]. In the vicinity of the caustic, the partial amplitude of HHG is described by Airy function, which forms a wide peak [30, 65, 66]. Caustic generated by electrons with first return at 21 fs < t_r < 23 fs leads to the formation of an intermediate 'plateau' in the HHG spectrum ending at ω ≈ 130 eV. High-order returns at 18 fs < t_r < 23 fs result in local maximum in HHG spectrum at ω ≈ 90 eV. Spectral caustics corresponding to the highest possible return energy lead to an increase in the spectral intensity in the cutoff region at 185 eV, and this increase is enhanced by the smooth growth of the differential PRCS in this frequency range.

4.1. Individual response of orbitals on two active shells

The TDDFT formalism allows one to separate the contributions of different KS orbitals to the total dipole acceleration, thereby providing a detailed study of their collective contribution to the enhancement of the HHG yield. In figure 4 we present the Fourier spectra of contributions a_j (t) for individual orbitals, of contributions ${\bar{a}}_{nl}\left(t\right)$ for subshells 4d, 4p, and 4s, and contribution for 4th shell determined by ${\hat{a}}_{n}\left(t\right)$ for n = 4.

Zoom In Zoom Out Reset image size

4.1.1. 5p₀ orbital

Let us consider the response of the outer-valence 5p₀ orbital, which is oriented in space along the electric field and interacts most actively with the linearly-polarized IR laser pulse. As seen from figure 4(a), the partial dipole acceleration spectrum of 5p₀ KS orbital contains a high-frequency flat plateau with the same cutoff as the full HHG spectrum. At low frequencies, ω ∼ 15–30 eV, the response of 5p₀ orbital contains a near-threshold enhancement associated with a sharp increasing of the transition matrix element from the continuum state to the 5p₀ state with decreasing frequency of emitted photons [11, 51, 53, 67, 68]. The response of the 5p₀ KS orbital is well described by the three-step scenario for the HHG [1, 8], i.e. the electron from the 5p₀ KS orbital is ionized by a laser pulse, then it is accelerated by a laser field in the continuum and, finally, electron recombines to the initial state with the emission of a photon with energy equal to the sum of gained kinetic energy and ionization potential. The only peculiarity is that response of 5p₀ KS orbital contains a narrow peak at ω = E_5p − E_4d ≈ 57.1 eV caused by the transition between 5p₀ and 4d₀ levels led by the electron–electron interaction. In the full HHG spectrum, this resonant peak is compensated by the corresponding resonance term in the response of 4d₀ KS orbital (see section 3). It should be emphasized that the spectrum of the dipole acceleration of 5p₀ KS orbital is similar to the HHG spectrum obtained with frozen electron–electron interaction approximation shown by the black line in figure 4(a). This explicitly shows that giant enhancement near 100 eV results from the strong interchannel coupling between different subshells.

4.1.2. Other orbitals from 5th and 4th shells

The tunneling rate in the infrared laser field of intensity 1.9 × 10¹⁴ W cm⁻² for orbitals others than 5p₀ is negligibly small. In the case of orbitals on the 4th shell and the 5s orbital, this is due to large binding energy (see table 1). In the case of 5p_±1 orbitals having the same binding energy as 5p₀, this is due to suppressing of the tunneling ionization for m ≠ 0 [69]. In addition, an important effect that reduces the probability of tunneling ionization for all KS orbitals, including 5p₀, is the polarization of the atom in an external field [20, 70, 71]. ⁸ As our calculations show, the final ionization probability of 5p_±1 orbitals is ≈0.6%; for other lower-lying orbitals, it is even less. Smallness in the tunneling rates suppresses the three-step scenario for all orbitals others than 5p₀. However, as our numerical calculation show (see figure 4(b)), 4d and 4p subshells also contribute to the total HHG yield. This non-three-step contribution can be explained by strong coupling of these subshells to the 5p₀ KS orbital and each other caused by electron–electron potential. Indeed, the time-dependent wave function of the 5p₀ KS orbital, ${\psi }_{5{p}_{0}}\left(\mathbf{r},t\right)\equiv {\psi }_{{j}_{5{p}_{0}}}\left(\mathbf{r},t\right)$, where ${j}_{5{p}_{0}}=26$, can presented as a sum of two terms [72–74]:

$$\begin{equation}{\psi }_{5{p}_{0}}\left(\mathbf{r},t\right)={\psi }_{\text{at}}\left(\mathbf{r},t\right)+{\psi }_{\text{res}}\left(\mathbf{r},t\right),\end{equation} \tag{ 4.4 }$$

where ψ_at(r, t) is so-called adiabatic part, which can be well approximated by the field-free wave function, ${\psi }_{\text{at}}\left(\mathbf{r},t\right)\approx {\psi }_{5{p}_{0}}^{\left(0\right)}\left(\mathbf{r}\right){\mathrm{e}}^{-\mathrm{i}{E}_{5{p}_{0}}t}$, and ψ_res(r, t) is the rescattering part of the wave function. Although this rescattering part is exponential small, it induces both the time-dependent moment given by the 5p₀ orbital, ${D}_{{j}_{5{p}_{0}}}\left(t\right)\approx \int {\psi }_{\text{at}}^{{\ast}}z{\psi }_{\text{res}}\enspace \mathrm{d}\mathbf{r}+\mathrm{c}.\mathrm{c}.$, and variation in the electron density ρ(r, t) from the corresponding field-free distribution,

$$\begin{equation}\rho \left(\mathbf{r},t\right)={\rho }_{0}\left(\mathbf{r}\right)+\delta \rho \left(\mathbf{r},t\right),\end{equation} \tag{ 4.5 }$$

$$\begin{equation}\delta \rho \left(\mathbf{r},t\right)\approx {\psi }_{\text{at}}^{{\ast}}{\psi }_{\text{res}}+\mathrm{c}.\mathrm{c}.\end{equation} \tag{ 4.6 }$$

The small changes in the density ρ(r, t) lead a variation of electron–electron interaction, δV_ee:

$$\begin{equation}{V}_{\text{ee}}\left[\rho \right]\approx {V}_{\text{ee}}\left[{\rho }_{0}\right]+\delta {V}_{\text{ee}},\hspace{25.0pt}\delta {V}_{\text{ee}}=\int \frac{\delta \rho \left({\mathbf{r}}^{\prime }\right)}{\vert \mathbf{r}-{\mathbf{r}}^{\prime }\vert }\mathrm{d}{\mathbf{r}}^{\prime }+{\dot {V}}_{\text{xc}}^{\left(0\right)}\delta \rho ,\end{equation} \tag{ 4.7 }$$

where ${\dot {V}}_{\text{xc}}^{\left(0\right)}\equiv \left(\partial {V}_{\text{xc}}/\partial \rho \right){\vert }_{\rho ={\rho }_{0}}$. The term δV_ee can be considered as a perturbation which acts as a broadband excitation source inducing the time-dependent moment D_j (t) of KS orbitals. Moreover, the rescattering part of 5p₀ KS orbital can affect on all orbitals including those with m ≠ 0. The account of δV_ee within perturbation theory leads to variation ψ_j (r, t) of the KS orbitals with index $j\ne {j}_{5{p}_{0}}$:

$$\begin{equation}{\psi }_{j}\left(\mathbf{r},t\right)\approx {\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right){\mathrm{e}}^{-\mathrm{i}{E}_{j}t}+\delta {\psi }_{j}\left(\mathbf{r},t\right),\end{equation} \tag{ 4.8 }$$

$$\begin{equation}\delta {\psi }_{j}=\sum\limits _{q}\delta {c}_{q}^{\left(j\right)}\left(t\right){\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right){\mathrm{e}}^{-\mathrm{i}{E}_{q}t},\end{equation} \tag{ 4.9 }$$

where the sum is carried out over all possible states q of the unperturbed Hamiltonian (2.1) forming a complete orthogonal system of wavefunctions ${\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right)$, and

$$\begin{align}\hfill \delta {c}_{q}^{\left(j\right)}\left(t\right)& =\delta {c}_{q}^{\left(j\right),+}\left(t\right)+\delta {c}_{q}^{\left(j\right),-}\left(t\right),\hfill \\ \hfill \delta {c}_{q}^{\left(j\right),+}\left(t\right)& =-\mathrm{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\Delta}{E}_{q,j}^{\left(+\right)}{t}^{\prime }}\iint \frac{{\left[{\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right){\psi }_{5{p}_{0}}^{\left(0\right)}\left({\mathbf{r}}^{\prime }\right)\right]}^{{\ast}}{\psi }_{\text{res}}\left({\mathbf{r}}^{\prime },{t}^{\prime }\right){\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)}{\vert \mathbf{r}-{\mathbf{r}}^{\prime }\vert }\mathrm{d}{\mathbf{r}}^{\prime }\enspace \mathrm{d}\mathbf{r}\hfill \\ \hfill & \quad -\mathrm{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\Delta}{E}_{q,j}^{\left(+\right)}{t}^{\prime }}\int {\left[{\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right){\psi }_{5{p}_{0}}^{\left(0\right)}\left(\mathbf{r}\right)\right]}^{{\ast}}{\dot {V}}_{\text{xc}}^{\left(0\right)}{\psi }_{\text{res}}\left(\mathbf{r},{t}^{\prime }\right){\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)\mathrm{d}\mathbf{r},\hfill \\ \hfill \delta {c}_{q}^{\left(j\right),-}\left(t\right)& =-\mathrm{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\Delta}{E}_{q,j}^{\left(-\right)}{t}^{\prime }}\iint \frac{{\left[{\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right){\psi }_{\text{res}}\left({\mathbf{r}}^{\prime },{t}^{\prime }\right)\right]}^{{\ast}}{\psi }_{5{p}_{0}}^{\left(0\right)}\left({\mathbf{r}}^{\prime }\right){\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)}{\vert \mathbf{r}-{\mathbf{r}}^{\prime }\vert }\mathrm{d}{\mathbf{r}}^{\prime }\enspace \mathrm{d}\mathbf{r}\hfill \\ \hfill & \quad -\mathrm{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\Delta}{E}_{q,j}^{\left(-\right)}{t}^{\prime }}\int {\left[{\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right){\psi }_{\text{res}}\left(\mathbf{r},{t}^{\prime }\right)\right]}^{{\ast}}{\dot {V}}_{\text{xc}}^{\left(0\right)}{\psi }_{5{p}_{0}}^{\left(0\right)}\left(\mathbf{r}\right){\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)\mathrm{d}\mathbf{r},\hfill \\ \hfill {\Delta}{E}_{q,j}^{\left({\pm}\right)}& ={E}_{q}-{E}_{j}{\pm}{E}_{5{p}_{0}}.\hfill \end{align} \tag{ 4.10 }$$

Using the variation in ψ_j , the time-dependent dipole moment of jth KS orbital (2.12) can be written as

$$\begin{equation}{D}_{j}\left(t\right)\approx \sum\limits _{q}{d}_{jq}\left(\delta {c}_{q}^{\left(j\right),+}\left(t\right)+\delta {c}_{q}^{\left(j\right),-}\left(t\right)\right){\mathrm{e}}^{\mathrm{i}\left({E}_{j}-{E}_{q}\right)t}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 4.11 }$$

where ${d}_{jq}=\int {\left[{\psi }_{j}^{\left(0\right)}\left(\mathbf{r}\right)\right]}^{{\ast}}z{\psi }_{q}^{\left(0\right)}\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}$ is the transition dipole matrix element. Among all possible intermediate states q, only those with the same quantum numbers m and differing by one orbital angular momentum l compared to the jth state contribute to the sum (4.11).

The high response of KS orbitals for j ≠ 5p₀ can be interpreted as a collective dynamics of the valence 5p₀ electron and a core electrons (see also references [14, 24, 31]). Indeed, in the first step, the valence electron (5p₀ electron) tunnels through the barrier and gets some extra kinetic energy E_k into the laser-perturbed continuum, which shares with a core electron on jth orbital. The correlation between these two electrons induces changing in electron–electron interaction and leads to the 'relaxation' of these two electrons into 5p₀ orbital and intermediate qth state with a hole in jth state. The subsequent transition from qth state to jth state returns the atom to the initial ground state with a harmonic emission of frequency $\omega ={E}_{k}-{E}_{5{p}_{0}}$ in agreement with the energy conservation law.

Excitation of oscillations of the 4th shell due to the recolliding photoelectron leads to the formation of a giant HHG yield in the frequency range ω = 60–130 eV, which at ω = 100 eV is more than one order of magnitude higher compared to the contribution of 5p₀ KS orbital. In this frequency range, the 4d subshell contains a giant resonance in PICS, so that 4d orbitals have large coupling matrix elements with the rescattering continuum wavepacket, 5p₀ state and intermediate f states in the continuum, and also high dipole matrix elements of transition f → 4d. Different orbitals from 4d subshell near the giant enhancement have the same high response (see figure 4(c) and reference [24]), which is due to the interaction of laser-perturbed 5p₀ KS orbital with all m orbitals and with the 'collectivization' of electrons in the 4d subshell [33, 40, 54]. However, as seen from figure 4(b), a large response in the region of giant HHG enhancement is also observed for deeper-lying 4p and 4s subshells. We note that it was previously accepted that only 4d subshell participates in HHG enhancement. Moreover, as can be seen from figures 4(c) and (d), the responses of individual KS orbitals in the 4th shell have a comparable response at ω > 100 eV. Coherent summation of responses from different subshells on 4th shell leads to a noticeable (more than two times) increase in HHG yield compared to the contribution of only 4d subshell in the frequency range ω ∼ 80–140 eV. Above the giant enhancement region, at ω > 170 eV, the response of the 4p subshell exceeds the response of the 4d subshell despite fewer electrons on the 4p subshell. The reason for the broadband high-frequency response of 4p and 4s orbitals is related, from the one hand, to the action of the V_ee perturbation induced by the rescattering wavepacket from 5p₀ orbital. From the other hand, as shown in section 3, 4p and 4s subshells can oscillate due to the coupling with the 4d orbital. In the case of HHG, this means that the recollision-induced perturbations δψ_j (4.9) in the wavefunctions of 4d KS orbitals induce a secondary perturbation in the electron–electron potential, which acts on all other KS orbitals, including 4p and 4s orbitals. Note that the effect of this secondary disturbance is also noticeable in the response of orbitals from 5th shell at ω = 70–130 eV (see figure 4(a)). In this frequency range signal amplification is observed for 5p_±1 orbitals. In general, the contribution of the 5p_±1 and 5s orbitals to HHG is significantly lower than that for the 4th shell and 5p₀ orbital because of the small matrix element of the transition to the high-energy continuum states. At the same time, induced time-dependent electron–electron potential δV_ee leads to a broadband response of the 5s and 5p_±1 KS orbitals on the outer shell for the frequency range ω ∼ 15–30 eV. A large low-frequency response of orbitals on the outer shell also exists in other noble-gas atoms, e.g. for Ar atom see references [20, 22]. The responses of orbitals on the outer shell are in phase and provide visible contribution (comparable with that for of 5p₀ KS orbital) to the total HHG yield for ω ∼ 15–30 eV.

To conclude, in this work, we have used the time-dependent density functional theory to study the HHG by Xe atom in an intense infrared pulse. In contrast to the previous numerical studies of this phenomenon, the dynamics of the whole 4th and 5th shells of the Xe atom is taken into account, leading to a good agreement with the experiment and the analytical parameterization of HHG yield in terms of the EWP and exact (experimental) PRCS to the 5p₀ state. It was demonstrated that as for argon atom [20], in the low-frequency range of generated spectral harmonics (ω ∼ 15–30 eV), a comparable contribution to HHG is given by all orbitals from the outer shell. The origin of such unusual collective contribution is the broadband component in electron–electron potential induced by the recolliding 5p₀ wave packet, which collisionally disturbs outer shell orbitals, thereby induces the time-dependent dynamics in these orbitals. In the giant enhancement region (ω ∼ 60–130 eV), the 5th shell's response is much smaller than that of the 4th shell, which is easily excited by the electron–electron potential due to the giant resonance in 4d transition to the continuum. Moreover, in addition to 4d subshell, deeper-lying 4p and 4s subshells also participate in the HHG, giving substantial contribution at higher frequencies. The reason for the substantial response of the 4p and 4s subshells at the giant resonance is the dynamic interaction between the orbitals on the 4th shell, associated with the close localization of their wave functions in space. The results of numerical simulations demonstrate that a similar response of 4p and 4s subshells is presented by interaction of Xe atom with an external linearly-polarized XUV field. Therefore, correlation interaction between all orbitals on the 4th shell can contribute to both the formation of the HHG spectra by infrared field and PICSs from 4d, 5p, and 5s subshells of Xe atom.

Numerical simulations were supported by the Russian Science Foundation (Grant No. 20-11-20289). This work was also supported by the Russian Foundation for Basic Research (Grant No. 20-32-70213), Grants Council of the President of the Russian Federation (Grant No. MK-1849.2020.2), and the Foundation for the Advancement of Theoretical Physics and Mathematics 'BASIS' (Grant No. 19-1-2-52-1).

All data that support the findings of this study are included within the article (and any supplementary files).

The HHG yield in an intense low-frequency linearly-polarized laser field can be factorized into two factors [30, 56, 57, 66, 75]. The first universal factor, $W\left(\mathcal{E}\right)$, known as EWP, describes tunneling and propagation of an active electron in a laser field along closed classical trajectories with returning energy $\mathcal{E}=\omega -{I}_{p}$, where ω is harmonic frequency. The second factor is field-free differential PRCS for the zero angle between returning electron momentum and polarization vector of emitted photon, ${\sigma }_{\text{rec}}\left(\mathcal{E},0\right)$,

$$\begin{equation}R\left(\omega \right)=\mathcal{W}\left(\mathcal{E}\right){\sigma }_{\text{rec}}\left(\mathcal{E},0\right).\end{equation} \tag{ A.1 }$$

Differential PRCS ${\sigma }_{\text{rec}}\left(\mathcal{E},\theta \right)$ is related to the differential PICS, σ_ion(ω, θ), by the principle of detailed balance [32]:

$$\begin{equation}{\sigma }_{\text{rec}}\left(\mathcal{E},\theta \right)={\sigma }_{\text{ion}}\left(\omega ,\theta \right)\frac{{\omega }^{2}}{2{c}^{2}\mathcal{E}},\end{equation} \tag{ A.2 }$$

where θ is an angle between photon polarization vector and electron momentum.

Although factorization (A.1) was theoretically justified within the single-active electron approximation (SAEA) [56, 66, 75], it provides a qualitative agreement with HHG experimental results [11, 12, 14–16]. For practical calculations, the form of EWP can be found numerically within SAEA calculation for a reference atom [56, 57], or either within analytical expressions [41, 76]. In this work, EWP was calculated analytically within adiabatic approximation for HHG [41]. Within adiabatic approximation, the EWP is determined by the set of real times {${t}_{i}^{\left(j\right)},{t}_{r}^{\left(j\right)}$}, corresponding to the moments of ionization (${t}_{i}^{\left(j\right)}$) and recombination (${t}_{r}^{\left(j\right)}$). These times obey a system of transcendental equations:

$$\begin{equation}{\boldsymbol{K}}_{j}^{\prime }=0,\hspace{25.0pt}\frac{{\boldsymbol{K}}_{j}^{2}}{2}=\omega -{I}_{p}-{\Delta}{\mathcal{E}}_{j},\end{equation} \tag{ A.3 }$$

where K '_j and K _j are instantaneous electron momenta at the moment of ionization and recombination

$$\begin{equation*}{\boldsymbol{K}}_{j}^{\prime }=\boldsymbol{K}\left({t}_{i}^{\left(j\right)};{t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right),\hspace{25.0pt}{\boldsymbol{K}}_{j}=\boldsymbol{K}\left({t}_{r}^{\left(j\right)};{t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right),\end{equation*}$$

and $\boldsymbol{K}\left(t;{t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)$ is instantaneous (at the moment t) momentum of an electron moving along a closed trajectory in the time interval $\left({t}_{i}^{\left(j\right)},{t}_{r}^{\left(j\right)}\right)$ in the laser field with vector potential $\mathbf{A}\left(t\right)=-\hat{\mathbf{z}}{\int }^{t}F\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }$,

$$\begin{equation}\boldsymbol{K}\left(t;{t}_{i}^{\left(j\right)},{t}_{r}^{\left(j\right)}\right)=\mathbf{A}\left(t\right)-\frac{1}{{t}_{r}^{\left(j\right)}-{t}_{i}^{\left(j\right)}}\underset{{t}_{i}^{\left(j\right)}}{\overset{{t}_{r}^{\left(j\right)}}{\int }}\mathbf{A}\left(\xi \right)\mathrm{d}\xi .\end{equation} \tag{ A.4 }$$

The term ${\Delta}{\mathcal{E}}_{j}$ is the quantum correction [72]

$$\begin{equation}{\Delta}{\mathcal{E}}_{j}=\frac{{I}_{p}}{{\Delta}{t}_{j}}\frac{{\boldsymbol{F}}_{j}^{\prime }\cdot {\boldsymbol{K}}_{j}}{{{\boldsymbol{F}}_{j}^{\prime }}^{2}},\quad \enspace \enspace \enspace \end{equation} \tag{ A.5 }$$

where ${\Delta}{t}_{j}={t}_{r}^{\left(j\right)}-{t}_{i}^{\left(j\right)}$ and ${\boldsymbol{F}}_{j}^{\prime }\equiv \boldsymbol{F}\left({t}_{i}^{\left(j\right)}\right)$. Detailed analysis of the system (A.3) for the case of linearly-polarized laser field can be found in [77].

The EWP $W\left(\mathcal{E}\right)$ is proportional to the square of a sum over the partial amplitudes, corresponding to closed electron trajectories,

$$\begin{equation}W\left(\mathcal{E}\right)=\frac{\pi }{2}\omega \sqrt{2\mathcal{E}}{\left\vert \sum\limits _{j}{a}_{j}^{\left(\mathrm{t}\mathrm{u}\mathrm{n}\right)}{a}_{j}^{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\right)}\right\vert }^{2},\end{equation} \tag{ A.6 }$$

where each partial amplitude can be factorized in terms of tunneling, ${a}_{j}^{\left(\mathrm{t}\mathrm{u}\mathrm{n}\right)}$, and propagation, ${a}_{j}^{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\right)}$, factors. The tunneling factor, ${a}_{j}^{\left(\mathrm{t}\mathrm{u}\mathrm{n}\right)}$, is given by

$$\begin{equation}{a}_{j}^{\left(\mathrm{t}\mathrm{u}\mathrm{n}\right)}={Q}_{j}\frac{{\mathcal{C}}_{1}}{\pi }\frac{{\mathrm{e}}^{-{\kappa }^{3}/\left(3{\mathcal{F}}_{j}\right)}}{\sqrt{2{\mathcal{F}}_{j}}},\end{equation} \tag{ A.7 }$$

where ${\mathcal{F}}_{j}=\sqrt{{{\boldsymbol{F}}^{\prime }}_{j}^{2}}$, $\kappa =\sqrt{2{I}_{p}}$ and ${\mathcal{C}}_{1}$ is dimensionless asymptotic coefficient of the initial bound state for the neutral atom,

$$\begin{equation*}{\left.{\psi }_{\kappa ,1,m}\right\vert }_{\kappa r\gg 1}\to {\mathcal{C}}_{1}\sqrt{\kappa }{r}^{-1}{\left(\kappa r\right)}^{1/\kappa }\enspace {\mathrm{e}}^{-\kappa r}{Y}_{1,m}\left(\theta ,\phi \right).\end{equation*}$$

The Coulomb factor, Q_j , in equation (A.7) for the case of linearly polarized laser field reduces to a 'static' factor [78],

$$\begin{equation}{Q}_{j}={Q}_{\text{stat}}^{\left(j\right)}={\left(\frac{2{\kappa }^{3}}{{\mathcal{F}}_{j}}\right)}^{1/\kappa }.\end{equation} \tag{ A.8 }$$

The propagation factor, ${a}_{j}^{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\right)}$, is given by

$$\begin{equation}{a}_{j}^{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\right)}\left({{\Omega}}_{h}\right)=\mathrm{i}\frac{{\mathrm{e}}^{-\mathrm{i}\mathcal{S}\left({t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)+\mathrm{i}\omega {t}_{r}^{\left(j\right)}}}{{\Delta}{t}_{j}^{3/2}\sqrt{{\boldsymbol{K}}_{j}\cdot {\dot {\boldsymbol{K}}}_{j}}},\end{equation} \tag{ A.9 }$$

where $\dot {{\boldsymbol{K}}_{j}}\equiv \partial \boldsymbol{K}\left({t}_{r}^{\left(j\right)};{t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)/\partial {t}_{r}^{\left(j\right)}$, and $\mathcal{S}\left({t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)$ is the classical action of an electron

$$\begin{equation}\mathcal{S}\left({t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)={I}_{p}{\Delta}{t}_{j}+\frac{1}{2}\underset{{t}_{i}^{\left(j\right)}}{\overset{{t}_{r}^{\left(j\right)}}{\int }}{\boldsymbol{K}}^{2}\left(\xi ;{t}_{r}^{\left(j\right)},{t}_{i}^{\left(j\right)}\right)\mathrm{d}\xi .\quad \end{equation} \tag{ A.10 }$$

The angular distribution of PICS follows from parametrization

$$\begin{equation}{\sigma }_{\text{ion}}\left(\omega ,\theta \right)=\frac{{\sigma }_{\text{ion}}\left(\omega \right)}{4\pi }\left[1+\beta \left(\omega \right){P}_{2}\left(\mathrm{cos}\theta \right)\right].\end{equation} \tag{ A.11 }$$

Here, σ_ion(ω) is the angle-integrated partial cross section for 5p subshell, β(ω) is the photoelectron asymmetry parameter, P₂(cos θ) = (3 cos²  θ − 1)/2 is a Legendre polynomial and θ is the angle between the photoelectron momentum direction and the photon polarization direction, ω is the photon energy. The values of σ_ion(ω) and β(ω) are taken from the results of experimental measurements in references [49, 62].