High-order harmonic generation (HHG) is generally believed to be the fundamental process underlying a broad range of attosecond applications, such as attosecond control and measurements,1–8 and is therefore one of the main focuses of current attosecond science. Increasing the brightness and cutoff energy,9–12 i.e., the maximum photon energy, of high-order harmonics has been of central importance for applications of HHG ever since its first demonstration in gases,2 since this enables higher resolution in photography of ultrafast electrons,13 higher accuracy in tracking the evolution of chemical reactions,14 and higher precision in steering electron wave packets in molecules.9,10,15

At present there are two parallel approaches to the generation of high-order harmonics, based on different theoretical principles. HHG in gases is widely used to produce coherent extreme-ultraviolet attosecond pulses.2,5,16 The mechanism involved here can be described in terms of a semiclassical “three-step” process16–19 in which high-order harmonic radiation is generated when an electron, ionized from its ground state, re-collides with its parent ion after it has oscillated for half a cycle in the intense driving laser field. The scenario for HHG in solids is quite different, involving several competing mechanisms, each of which provides a channel for HHG.20–25 In addition to interband transitions, which are reminiscent of the atomic semiclassical mechanism, intraband Bloch oscillations of electrons and holes derived from lattice structure are also important sources of high-order harmonics in solids.26–28 For the most recent progress in this direction, the reader should consult Ref. 29 and the description therein of the Wannier quasiclassical model, which illustrates the importance of intersite interaction. Emitted high-order harmonic radiation is the result of interfering superposition of all these channels.

From the perspective of applications of HHG, the “three-step” mechanism in gases might be considered preferable, since it provides a single channel of HHG and therefore much easier to control.17,19 However, as far as the intensity of emitted radiation is concerned, the high density of atoms in solids and their regular lattice configuration are of particular significance.20,30,31 It has been demonstrated experimentally by Popmintchev et al.11 that a fine balance between a high intensity of emitted radiation and the atomic mechanism of HHG can be attained in gases at high pressure (up to 30 atm for noble gases), corresponding to an average distance d₀ between two atoms of about 11 Å. At this distance, there is negligible overlap of the ground-state wave functions of different atoms, which effectively suppresses the HHG channel associated with the hole dynamics typical of solids. However, experiments have also revealed that a new channel of HHG emerges, associated with re-collisions of electrons with nearby ions.

How the HHG mechanism is transformed from the atomic “three-step” one17,19 to the multichannel one32–34 in solids is therefore an important question, the answer to which may help provide not only further insights into HHG itself, but also detailed theoretical guidance for the design of HHG devices. In addition to various simplified theoretical models of the transition of the HHG mechanism,18,33,35–38 there are also some first-principles approaches that have been applied successfully to this problem. One of these is the time-dependent density functional theory (TDDFT) method,39–41 which affords a quantum-mechanical description of the time evolution of electrons in the system based on the time-dependent Kohn–Sham equations. Previous work on semiconductors and insulators has shown that the TDDFT method performs satisfactorily in its predictions of HHG in various solids.32,42–47 However, it is not trivial to apply the same approach to the description of HHG in a system of gaseous density, because of the tremendous computational cost then incurred, since this cost scales roughly as d03, with d₀ being the distance between atoms. However, with the most recent advances in computational techniques and a much improved implementation of the TDDFT method,48–51 it is now possible to get reliable predictions for low-density systems such as noble gases at 0.06 atm, corresponding to d₀ ≈ 85 Å (160 bohrs). This paves the way for the creation of a physical picture of the transition of the HHG mechanism.

In this work, with the TDDFT method39–41 as our tool, we present a systematic numerical study of the transition of HHG in neon (Ne) systems. The density of the system varies from that of an isolated Ne atom to that of solid Ne, namely, 3.55 × 10²² cm⁻³, corresponding to an average distance d₀ = 3.04 Å (5.75 bohrs). We show that in an extended system, ionized electrons wander about in the system until they meet a nearby ion and collide with it. There is no maximum distance limit as in the atomic HHG mechanism, and the maximum kinetic energy of the recombined electron is then determined to be 3.2U_p,17,52 with U_p = I₀/(4ω²) being the ponderomotive energy delivered by a laser of frequency ω and intensity I₀. Note that throughout this article, the formulas are all expressed in atomic units. The edge of the main feature of the high-order harmonic spectrum in this “wandering electron” picture is then given by E_edge = I_p + 8U_p, where I_p is the ionization energy of Ne. The factor of 8 comes from the maximum kinetic energy of an ionized electron in the driving laser field. Beyond the atomic limit of high-order harmonic spectra, a multiplatform feature is observed, corresponding to re-collisions of ionized electrons with nearby ions. In addition, it is shown that a Ne simple cubic lattice of size 63.5 Å (120 bohrs) provides a selection condition for the direction of polarization of high-order harmonics beyond the atomic limit, which may be further used to manipulate the emitted radiation.

The remainder of this article is organized as follows. After reviewing the theoretical method in Sec. II, we present our main results in Sec. III, where the transition of the HHG mechanism and the energy edge in a low-density gas system are discussed in detail. A short summary is provided in Sec. IV.

TDDFT describes the evolution of Kohn–Sham orbitals of a multi-electron system using the time-dependent Kohn–Sham equations,i∂∂tψα(r,t)=−∇22+Vext(r,t)+VH[n(r,t)]+Vxc[n(r,t)]ψα(r,t),(1)following the arguments of Runge and Gross.39–41 Here, α is the index of a Kohn–Sham orbital, ψ_α(r, t) is the wave function, V_ext is the external potential (including the ionic potential, the trapping potential of the atoms, and the applied laser field), V_H is the Hartree part of the Coulomb interaction between electrons, and V_xc is the time-dependent exchange-correlation potential. The electron density is n(r, t) = ∑_α|ψ_α(r, t)|², and the electronic current density j(r, t) is calculated as32,53j(r,t)=i2∑αψα*(r,t)∇ψα(r,t)−ψα(r,t)∇ψα*(r,t).(2)The high-order harmonic spectrum is obtained as the Fourier transform of the rate of change of the laser-driven electronic current:32HHG(ω)=FT∂∂t∫Ωd3rj(r,t)2,(3)where Ω represents the volume of the system and FT denotes the Fourier transform operation. Note that for the calculation of the time-dependent high-order harmonic spectrum, a small Gaussian time window of 0.13 fs width will be used in Eq. (3) before Fourier transformation. When attosecond pulses [e.g., an isolated attosecond pulse (IAP) or a pair of attosecond pulses] can be constructed from high-order harmonic spectra, they are calculated as coherent superpositions of consecutive harmonics in a given frequency window from ω_l to ω_h as follows:9I(t)=∑ωlωheiωtFT∂∂t∫Ωd3rj(r,t′)2.(4)

In our calculations, Ne atoms are placed at lattice points of a three-dimensional simple cubic optical lattice for simplicity of calculation. The optical lattice can be generated by the potential of three pairs of identical laser beams,54 as illustrated schematically in Fig. 1(a), which form standing waves in three orthogonal directions. The resulting optical potential can be written asV(x,y,z)=V0[sin2(kLx)+sin2(kLy)+sin2(kLz)],(5)which provides the symmetry of a simple cubic lattice. A top view of the potential is shown in Fig. 1(b). Here k_L is the wave vector of the periodic potential, determined by the wavelength λ_L of the laser beam as k_L = 2π/λ_L. The lattice constant is given by d = λ_L/2. V₀ is the potential depth of the standing waves, sometimes expressed using the recoil energy Er=kL2/(2m) as the unit, with m being the mass of the Ne atom, as a direct measure of its ability to confine atoms. In our calculations, E_r is of the order of 10⁻⁴ eV, and V₀ is set to be about 100E_r and thus has negligible influence on the energetic behavior of electrons. It should be noted that with current techniques, it is only possible to generate optical lattices with lattice constants of several hundreds of nanometers or longer. Optical lattices with a shorter lattice constant are still theoretical constructs. However, they will be considered as a convenient theoretical tool in this work, mainly to take account of density effects.

All calculations are performed using the OCTOPUS code,48–51 which is a real-space implementation of the TDDFT method. A primitive unit cell with periodic boundary conditions is used, except for the case of an isolated Ne atom, where a spherical calculation region of 100 bohrs (52.9 Å) radius is used with absorbing boundary conditions. The width of the absorbing region is 25 bohrs (13.2 Å), inside which there is a sine-square-shaped imaginary absorbing potential with a scaling factor −2.0 hartree (−54.4 eV), following De Giovannini et al.55 The resolution of the real space is 0.4 bohr (0.21 Å), and the energy cutoff of the kinetic energy is set to be 30.84 hartree (839.2 eV). The k-point mesh used to sample the Brillouin zone56 varies from 2 × 2 × 2 for d = 160 bohrs (84.7 Å) to 24 × 24 × 24 for d = 5.75 bohrs (3 Å), increasing with decreasing lattice constant. The time step to advance the evolution of the electronic wave functions is 0.48 as, and the exchange-correlation interaction is described by the adiabatic local density approximation (ALDA), with the parameterization of Perdew and Zunger.57 The Ne pseudopotential has a norm-conserving Troullier–Martins58 form with eight electrons included explicitly and a cutoff radius of 2.6 bohrs (1.38 Å). Before the system is put into time-dependent evolution, the initial structure and ground-state electronic structure are first generated from a time-independent density functional theory (DFT) calculation.59,60

The driving laser pulse is described by a vector potentialA(t)=8πI0/cωS(t)cos(ωt+ϕ0)ê,(6)where I₀ is the peak intensity, c is the speed of light, ω is the fundamental frequency of the driving laser, S(t) is the envelope function, ê is the polarization vector, and ϕ₀ is the carrier envelope phase. In the calculations, the fundamental wavelength of the driving laser is 8000 Å, corresponding to ω = 1.55 eV, and the peak intensity is between 50 and 700 TW/cm². In the calculation of high-order harmonic spectra, the pulse has a duration T = 16 fs, corresponding to six oscillating cycles, and a sine-square envelope S(t) = sin²(πt/T), as displayed in Fig. 2(a). In the calculation of the IAP, the envelope is set as S(t) = sin^3.7(πt/T), with a pulse duration of 13.33 fs, i.e., five oscillating cycles. The carrier envelope phase ϕ₀ is set to be 3π/2 for the generation of the IAP; otherwise, it is set to be zero. Since the wavelength of the driving laser is much longer than the lattice constant d, which varies from 5.75 bohrs (3 Å) to 160 bohrs (84.7 Å), the field felt by the electron in the primitive cell is approximated as a time-dependent uniform field.

We first show the general features of HHG in a Ne system at a density comparable to gaseous density. For simplicity of calculation, the Ne atom is assumed to be placed in a simple cubic optical lattice. In Fig. 2(b), a typical high-order harmonic spectrum of Ne is displayed (on a logarithmic scale), where the intensity of the laser is 500 TW/cm², the driving field is along the [100] direction of the simple cubic lattice, and the lattice constant d = 120 bohrs (63.5 Å), corresponding to a number density of 3.91 × 10¹⁸ cm⁻³. The pulse duration is 16 fs, consisting of six oscillating cycles of the driving laser, as displayed in the inset of Fig. 2(a). A time-dependent high-order harmonic spectrum, calculated as the Fourier transform in a small time window of 0.13 fs, is also displayed in Fig. 2(d) to show the timing of HHG.

For comparison, the high-order harmonic spectrum of an isolated Ne atom is also displayed in Fig. 2(a). The atomic energy cutoff of high-order harmonics given by Refs. 17 and 52E_cut = I_p + 3.2U_p is shown by the vertical blue dashed line. It can be seen from Fig. 2(a) that at E_cut, which is about 110 eV or the 70th order of harmonics for an Ne atom, the spectrum displays a sharp edge, beyond which the intensity decreases to about 10⁻⁴ of its maximum value. The time-dependent high-order harmonic spectrum in Fig. 2(c) shows that higher-order harmonics are generated in every half cycle of the driving laser field. Both features agree well with the prediction of the semiclassical “three-step” model, showing that the TDDFT calculation indeed captures the main features of HHG of isolated atoms.

By contrast, the high-order harmonic spectrum of an extended system has a typical multiplatform structure, as displayed in Fig. 2(b). Below the atomic energy cutoff, the spectrum is very similar to that of an isolated atom, showing the persistence of the semiclassical mechanism of HHG, i.e., the recombination of ionized electrons with their parent ions. Beyond the energy limit, two extra platforms are observed, corresponding to the re-collision of electrons with their nearest-neighbor and next-nearest-neighbor ions. The platform region ends up at a new energy limit of about 250 eV, corresponding to the 160th order of harmonics. As we shall show in the following discussion, this energy limit is determined by the maximum energy delivered by the driving laser to an ionized electron.

For isolated atoms, the energy cutoff E_cut of the high-order harmonic spectrum is mainly determined by the maximum distance that an electron can reach before re-collision with its parent ion.17,52 In an extended system, however, one can expect that this distance limit will no longer apply. Instead, in a semiclassical picture, an ionized electron can wander around the extended system until it meets an ion and recombines with it. Although the position of the electron is not known exactly during its wandering, its maximum kinetic energy is determined by the acceleration and deceleration process of the driving field and turns out to be 8U_p at the moment when the electron is accelerated for a half cycle. Thus, the maximum energy of the high-order harmonic spectrum is given by the maximum kinetic energy plus the ionization energy, i.e., Ecutnew=Ip+8Up.

The wandering electron picture suggests that the HHG process is no longer synchronized with the driving cycles, and this is indeed observed in the time-dependent high-order harmonic spectrum displayed in Fig. 2(d). Compared with the time-dependent high-order harmonic spectrum of an isolated Ne atom in Fig. 2(c), it can be seen that the high-order harmonics can still be generated in the fifth cycle. Since the laser field in the fifth cycle is not strong enough to ionize additional electrons, as illustrated in the case of the isolated Ne atom, these high harmonics are actually generated by wandering electrons ionized in previous driving cycles.

To confirm that the platforms beyond the atomic energy cutoff in Fig. 2(b) originate from the collision of an electron with its nearby ions, energy distributions of an electron at the times at which it arrives at the nearest-neighbor ion and the next-nearest-neighbor ion are plotted in Fig. 3. The results are calculated using the semiclassical model, which assumes that an electron is accelerated only by the driving field without an initial velocity after it was ionized at a certain moment in the driving cycle. The ionization probability is assumed to be proportional to the intensity of the driving field, which is set to be monochromatic with an intensity of 500 TW/cm² at λ_L = 8000 Å. The lattice constant used in the calculation is d = 120 bohrs (63.5 Å). Figure 3 shows that the energy distribution of an electron arriving at the nearest-neighbor ion (dark cyan curve) has a peak at around 130 eV, and goes down to zero at the maximum energy 250 eV. The energy distribution of an electron at the next-nearest-neighbor ion (violet curve) has a peak at a higher energy around 200 eV. The total distribution (orange curve) agrees qualitatively with the TDDFT result displayed in Fig. 2(b).

To understand how the transition in the HHG mechanism occurs with increasing ionic density, a series of calculations are carried out at different ionic densities ranging from nearly zero, i.e., the isolated atom case, to the density of solid Ne, where d is about 5.75 bohrs (3 Å). Figure 4 shows the calculated high-order harmonic spectra for several ionic densities. The density of the system decreases from Figs. 4(a)–4(d). It can be seen that below a density of 1.65 × 10¹⁸ cm⁻³, where d = 160 bohrs (84.7 Å), the high-order harmonic spectra are atomic-like. The wandering of electrons is not significant. At higher densities, the high-order harmonic spectra show typical multiplatform features beyond the atomic energy cutoff. All of these features disappear at the maximum energy around 8U_p + I_p. It is interesting to notice that, in this density regime, the number of platforms increases with increasing density, showing that the electrons can recombine with ions that are not close to their parent ions, which provides further support for the picture of wandering electrons. When the ionic density reaches 3.55 × 10²² cm⁻³, i.e., the density of solid Ne, the high-order harmonic spectrum still displays some of the features of wandering electrons. In particular, the maximum energy limit of 8U_p + I_p can still be observed, but the multiplatform features are smoothed out.

The transition from an atomic-like spectrum of an isolated atom to the multiplatform spectrum of a low-density extended structure in Fig. 4 is due mainly to the finite duration of the driving laser pulse. Figure 5 shows the energy distributions of an electron ionized at different moments in the first oscillating cycle of the driving laser. The distributions are calculated at the ends of the first, second, and third cycles, respectively, at each of which two peaks are observed, corresponding to the energy distribution of electrons ionized in the first and second half cycles, respectively. The three distributions show that the maximum energy of the electron recorded at the end of an oscillating cycle is less than 70 eV, which is below the atomic energy limit of 110 eV. This suggests that the energy of an electron will not exceed the atomic energy limit given by Refs. 17 and 52I_p + 3.2U_p if the pulse finishes with a complete driving cycle. Therefore, for a very long distance between two nearby atoms, the ionized electron is not able to reach its nearest-neighbor ions before the end of the pulse duration, and the energy it carries will not exceed the atomic energy limit, as revealed in Fig. 5. At the given laser intensity of I_L = 500 TW/cm², the transition takes place at a density of about 1.65 × 10¹⁸ cm⁻³, corresponding to d = 160 bohrs (84.7 Å).

The transition of the high-order harmonic spectrum to one of a solid type at further increased density probably results from the formation of a solid-like band structure,61 where significant band distortion occurs at the edges of Brillouin zones owing to strong interactions between atoms at short distances. When this condition holds, the HHG mechanism is dominated by a combination of Bloch oscillations and electronic transitions between bands.34,62 However, the transition of the band structure is not significant even in the solid Ne case. The band structure of the simple cubic structure at the density of solid Ne is displayed in Fig. 6, where it can be seen that the band structure is quite close to that of a free electron, shown by dotted curves, as a result of weak van der Waals interactions between Ne atoms. Therefore, it is not surprising that the high-order harmonic spectrum at this density still has the same maximum energy limit of about 250 eV.

One can also calculate the dependence of the high-order harmonic energy edge on the intensity of the driving laser. Figure 7 shows the calculated intensity dependence at an ionic density of 3.91 × 10¹⁸ cm⁻³, corresponding to d = 120 bohrs (63.5 Å). It can be seen that there is a minimum intensity threshold at around 50 TW/cm² for the high-order harmonics being observed. Furthermore, there is a threshold intensity of about 250 TW/cm², below which the maximum energy of high-order harmonics satisfies the atomic limit, as a result of finite pulse duration, i.e., the edge of the main feature is I_p + 3.2U_p. When the intensity is above the threshold, the energy edge closely follows the straight line corresponding to 8U_p + I_p.

The coherence of the high-order harmonics is retained well in the low-density extended system. A good indication of this is the generation of attosecond pulses. The inset of Fig. 2(b) shows the generation of a pair of attosecond pulses, which are calculated in a simple cubic periodic system with d = 120 bohrs, a laser intensity I_L = 500 TW/cm², and a sine-square-shaped driving laser pulse containing six cycles, as described in Sec. II. The carrier envelope phase ϕ₀ is set to be zero, and the energy window is between 158.2 and 300 eV [the interval between the two vertical red dashed lines in Fig. 2(b)]. From the time-dependent high-order harmonic spectrum in Fig. 2(c), it can be seen that the harmonics in the given energy window start from the third cycle of the pulse and last for about three cycles, which suggests that the attosecond pulses result from coherent superposition of high-order harmonics with interference cancellation. The width of each pulse is about 60 as, and the separation between the two pulses is about half of the cycle, i.e., 1.3 fs. An IAP can also be generated from high-order harmonics with a shorter driving pulse and an appropriately adjusted carrier envelope phase. Figure 8 displays an IAP calculated with a femtosecond driving pulse, for which the envelope function S(t) = sin^3.7(πt/T), with T = 13.33 fs, i.e., five oscillating cycles. The energy window is between 140 and 300 eV. In this way, it is possible to generate an IAP with a width of about 57 as.

The wandering electron picture and the multiplatform feature of the high-order harmonic spectrum are associated in particular with the effect of the density of the system. However, an extra dimension for the manipulation of HHG may be afforded by the presence of a periodic structure. Figure 9 illustrates how the direction of polarization can be selected in a simple cubic lattice. This selection takes advantage of the fact that the transition from an isolated atomic high-order harmonic spectrum to the multiplatform spectrum occurs at a lattice constant d ≈ 160 bohrs (84.6 Å) for a laser intensity I_L = 500 TW/cm². When the lattice constant d is set to be 120 bohrs (63.5 Å), the lattice can exert a selection effect on the direction of polarization of the driving laser. As can be seen in Fig. 9, only when the polarization is parallel to the [100] direction can the high-order harmonics beyond the atomic energy limit be generated. In the other two diagonal directions, i.e., along [110] and [111], the multiplatform feature is suppressed because of the increased distance between Ne atoms.

Using real-time time-dependent density functional theory calculations, we have systematically investigated the transition of the HHG mechanism in an isolated Ne atom to that in an extended Ne system of solid density. For the generation of high-order harmonics in such systems, we have presented a wandering electron picture that accounts for the multiplatform feature of high-order harmonic spectra and predicts the energy edge to be 8U_p + I_p. We have also shown that the high-order harmonics generated in this way may be used to generate an IAP of pulse width about 57 as. In addition, we have found that periodic lattices exhibit an extra selection effect depending on the direction of polarization of the driving laser, which may provide a new way to manipulate high-order harmonics. These interesting observations may help to provide not only a more complete picture of the HHG mechanism, but also new theoretical insights to aid in the design of future HHG devices.

Acknowledgment. This work is financially supported by the National Key Research and Development Program of China (Grant No. 2017YFA0404200) and the National Natural Science Foundation of China (Grant Nos. 12174034 and 91752202). Some of the calculations were carried out at the National Supercomputer Center in Gangzhou.