Determining the optimum thickness for high harmonic generation from nanoscale thin films: an ab initio computational study
DDetermining the optimum thickness for high harmonic generation from nanoscale thinfilms: an ab initio computational study
Shunsuke Yamada and Kazuhiro Yabana Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan (Dated: March 1, 2021)We theoretically investigate high harmonic generation (HHG) from silicon thin films with thick-nesses from a few atomic layers to a few hundreds of nanometers, to determine the most efficientthickness for producing intense HHG in the reflected and transmitted pulses. For this purpose, weemploy a few theoretical and computational methods. The most sophisticated method is the abinitio time-dependent density functional theory coupled with the Maxwell equations in a commonspatial resolution. This enables us to explore such effects as the surface electronic structure andlight propagation, as well as electronic motion in the energy band in a unified manner. We alsoutilize a multiscale method that is applicable to thicker films. Two-dimensional approximation isintroduced to obtain an intuitive understanding of the thickness dependence of HHG. From these abinitio calculations, we find that the HHG signals are the strongest in films with thicknesses of 2–15nm, which is determined by the bulk conductivity of silicon. We also find that the HHG signalsin the reflected and transmitted pulses are identical in such thin films. In films whose thicknessesare comparable to the wavelength in the medium, the intensity of HHG signals in the reflected(transmitted) pulse is found to correlate with the magnitude of the electric field at the front (back)surface of the thin film.
I. INTRODUCTION
Following progress in high harmonic generation (HHG)in atoms and molecules that enabled the productionof attosecond pulses , HHG in bulk crystals has at-tracted great interest during the last decade . Ex-perimental studies have achieved HHG in thin films ofvarious materials and laser pulses . These are ex-tended to monatomic two-dimensional (2D) layers, suchas graphene , and to metasurfaces, periodic 2D struc-tures composed of nanoscale objects . HHG from bulkcrystals is expected to be very intense compared withthose from atoms in the gas phase because of the highatomic density. This indicates that the HHG from solidsis favorable in applications to develop, for example, com-pact devices of XUV light sources.Theoretical studies in the past mostly focused on elec-tronic motion induced by a pulsed electric field preparedin advance . For sufficiently thick films, it is legiti-mate to consider electronic motion in an infinitely pe-riodic crystalline system in three dimensions induced bya spatially uniform electric field. Calculations employingtheories of varying complexity, such as one-dimensionalmodel, , time-dependent Schr¨odinger equation ,density matrix models , Floquet theory and abinitio descriptions such as the time-dependent densityfunctional theory (TDDFT) , have been developed.Theories to describe HHG from monatomic layers andbulk surfaces have also been developed . Using thesetheories, electronic motion in the wavenumber ( k ) spacebased on the energy band picture has been investigatedand classified as inter- and intra-band motions . It hasbeen clarified that there are rich manifestations of thecrystalline structure of HHG through the band struc-ture, selection rules, and dependence on the polarizationdirection . To investigate HHG in bulk materials, the propaga-tion effect is also important. The light-propagation ef-fect on HHG has been considered already in the case ofatomic gases . However, the effect should be significantprimarily in bulk solids, owing to their high atomic den-sity. The propagation effect in HHG from thin films hasbeen investigated recently. Experimentally, HHG in re-flected and transmitted pulses has been measured andcompared for films of several different thicknesses .Theoretically, multiscale calculations coupling a coarse-graining electronic motion with light propagation havebeen conducted . It has been shown that thelight propagation works to produce a clear HHG spec-trum, even when the spectrum is unclear in the unit-cell calculation . It has also been shown that the HHGin the transmitted pulse decreases by several orders ofmagnitude when emitted from thin films of micro-meterthickness .Although there has been progress on the propagationeffect, as described above, we consider that even the fol-lowing basic questions have not yet been answered: (1)At what film thickness is HHG emitted with the maxi-mum intensity? (2) What are the differences and similar-ities between HHGs in reflected and transmitted pulses?The principal purpose of this paper is to answer thesequestions unambiguously.We theoretically describe HHG in Si thin films of var-ious thicknesses, from a few atomic layers to a few hun-dreds of nanometers, that are comparable to the wave-length of the incident pulse in the medium. We utilize afew theoretical and computational frameworks based on ab initio TDDFT coupled with the Maxwell equa-tions. In the most sophisticated method, we simulta-neously solve the time-dependent Kohn-Sham (TDKS)equation for the electronic motion in TDDFT and theMaxwell equations for electromagnetic fields with a com- a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b mon spatial resolution. We call this the single-scaleMaxwell-TDDFT method . The method is applica-ble to films with thicknesses of less than a few tens ofnanometers because of its high computational cost. Forthicker films, we utilize a multiscale Maxwell-TDDFTmethod, in which a microscopic TDKS equation is cou-pled with the macroscopic Maxwell equations with acoarse-graining approximation . A 2D approximation is also introduced to obtain an intuitive understandingof the thickness dependence of HHG. Calculations usingthese methods provide a unified understanding on theHHG in nanoscale thin films.This paper is organized as follows: Sec. II describesthe theoretical and computational methods based onTDDFT. In Sec. III, the calculation results are presented.The thickness dependence of HHG, as well as the relationbetween HHGs in reflected and transmitted pulses, arediscussed. Finally, a conclusion is presented in Sec. IV. II. THEORETICAL METHODSA. Problem setup and method summary
We consider irradiation of a free-standing thin film ofthickness d in a vacuum by an ultrashort light pulse of alinearly polarized plane wave at normal incidence. As aprototypical material, we selected Si thin films of variousthickness, from a few atomic layers to a few hundreds ofnanometers. We assume that the thin film is infinitelyperiodic and macroscopically isotropic in 2D. The prop-agation and polarization directions of the incident lightare taken to be along the z and x axes, respectively. Weassume that the reflected and transmitted pulses containonly a component in the x direction, owing to symmetry.We express the asymptotic form of the vector potentialin the following way using the x component of the vectorpotential A ( z, t ): A ( z, t ) = (cid:40) A (i) ( t − z/c ) + A (r) ( t + z/c ) ( z → −∞ ) ,A (t) ( t − z/c ) ( z → + ∞ ) , (1)where A (i) ( t ), A (r) ( t ), and A (t) ( t ) are the incident, re-flected, and transmitted fields, respectively. The x com-ponent of the electric field, E ( t ), is related to the vec-tor potential by E ( t ) = − (1 /c ) dA ( t ) /dt . The reflectedand transmitted electric fields, E (r) ( t ) and E (t) ( t ), con-tain the high-order harmonics generated by the nonlinearlight-matter interaction in the film.To calculate HHG spectrum from a thin film theoreti-cally, it has been commonly considered electronic motionin a unit cell of crystalline solids induced by an appliedelectric field whose time profile is prepared in advance.There are, however, a number of effects that should beconsidered further. We investigate the following threeeffects in this paper. First is an effect of the surface elec-tronic structure, that is, a change of the electron bandstructure at the surfaces from that of the bulk system. The other two effects are related to the difference betweenthe incident electric field E (i) ( t ) and the electric field thatactually acts on electrons in the thin film. We considertwo effects that cause this difference. One is the effectof the light propagation. A strong incident laser pulse ismodulated during the propagation due to nonlinear light-matter interactions, and HHG pulses that are generatedinside the thin film also suffer from strong modulationand absorption during propagation before they exit themedium at the surfaces. The other is the effect of theconductivity of the thin film. Even when the thin film isso thin that the electric field inside the thin film can beregarded as uniform at the macroscopic scale, the electricfield inside the film is different from that of the incidentpulse due to the conductivity. Specifically, the currentthat flows in the film causes reflection from the film, andthe electric field inside the medium is equal to the sum ofthe incident and reflected electric fields, not the incidentfield alone.For a theoretical and computational description ofHHG from a thin film, we utilize four ab initio approachesbased on TDDFT , which are summarized in Table Iand explained in the following subsections. The simplestapproach is the single unit-cell method, in which we solvethe TDKS equation for electronic motion in a unit cellof the crystalline solid. In this description, the electricfield of the incident pulse is used as an applied field. Themost extensive and sophisticated approach is the single-scale Maxwell-TDDFT method , in which we solve theMaxwell equations for the electromagnetic fields and theTDKS equation for the electronic motion simultaneouslyusing a common spatial grid. Here, all three effects men-tioned above are included. In the multiscale Maxwell-TDDFT method , macroscopic light propagation andmicroscopic electronic motion are coupled using a coarse-graining approximation. In the practical calculation, theMaxwell equations and TDKS equation are solved simul-taneously using different spatial grids. While the effectof the surface electronic structure is not included, thismethod can describe HHG from thick films that cannotbe accessed by the single-scale Maxwell-TDDFT method.Finally, the 2D approximation is introduced as an ap-proximation to the multiscale Maxwell-TDDFT method,treating the thin film as an infinitely thin in macroscopicscale. This method will be useful to understand the thick-ness dependence of the HHG. B. Single unit-cell method
First, we consider a method describing the electronicmotion in a unit cell of a crystalline solid under aspatially uniform electric field . This method de-scribes electronic motion based on the energy-band pic-ture, and has been utilized extensively to describe phe-nomena related to nonlinear and ultrafast dynamics ofelectrons such as nonlinear susceptibilities , nonlinearenergy transfer , attosecond dynamics , and HHG TABLE I. Methods that will be used to describe HHG and their effectiveness to take account of various effects.Theory Single-scale Maxwell-TDDFT Multiscale Maxwell-TDDFT 2D approximation Single unit-cellTreatment of EM fields Microscopic Macroscopic 2D macroscopic NoneBulk electronic structure Yes Yes Yes YesSurface electronic structure Yes No No NoLight propagation Yes Yes No NoConduction Yes Yes Yes No from various solids . Within the dipole approxima-tion, the electron dynamics are described using the Blochorbital of the bulk system, u bulk n k ( r , t ), where the wavevector k runs over the three-dimensional (3D) Brillouinzone of the unit cell. The TDKS equation for the Blochorbitals u bulk n k ( r , t ) is written as i (cid:126) ∂u bulk n k ( r , t ) ∂t = (cid:26) m (cid:16) − i (cid:126) ∇ + (cid:126) k + ec ˆ x A (i) ( t ) (cid:17) − eφ ( r , t ) + δ ˆ V ion + V xc ( r , t ) (cid:27) u bulk n k ( r , t ) . (2)Here, we use the vector potential of the incident pulse, A (i) ( t ), as the applied electric field. The electricscalar potential φ ( r , t ) includes the Hartree potentialfrom the electrons and the local part of the ionic po-tential. δ ˆ V ion and V xc ( r , t ) are the nonlocal part ofthe ionic pseudopotential and the exchange-correlationpotential , respectively. All of the calculations withinthis paper are performed within the adiabatic local-density approximation . We ignore the exchange-correlation effect for the vector potential, A xc ( t ), forsimplicity .The electric current density averaged over the unit cellis calculated as follows: J ( t ) = − em (cid:90) Ω d r Ω occ (cid:88) n k u bulk ∗ n k ( r , t ) × (cid:16) − i (cid:126) ∇ + (cid:126) k + ec ˆ x A (i) ( t ) (cid:17) u bulk n k ( r , t ) + δ J ( t ) , (3)where Ω is the volume of the unit cell and the index n runs over the occupied band in the ground state. δ J ( t )is the current density from the nonlocal part of thepseudopotential . δ J ( t ) = − em (cid:90) d r Ω occ (cid:88) n k u bulk ∗ n k ( r , t ) e − i ( k +( e/c )ˆ x A (i) ( t )) r × [ r , δ ˆ V ion ] i (cid:126) e i ( k +( e/c )ˆ x A (i) ( t )) r u bulk n k ( r , t ) . (4)In practice, the current contains only the x component,parallel to the polarization direction of the incident pulse.We later discuss how to relate this current density withthe HHG spectrum from a thin film of thickness d . C. Single-scale Maxwell-TDDFT method
We next consider two theoretical methods that arecapable of describing the electronic motion and lightpropagation simultaneously. The first is the single-scaleMaxwell-TDDFT method, in which the electromagneticfields are treated microscopically. In the next subsec-tion, we will present the second method, the multiscaleMaxwell-TDDFT method, in which the electromagneticfields are treated macroscopically.The single-scale Maxwell-TDDFT method can fullyaccount for the effects of the light propagation, conduc-tion, and surface electronic structure. In this method,the Maxwell equations for electromagnetic fields and theTDKS equation for electronic motion are solved simulta-neously in the time domain using a common spatial grid.In the present case, we consider an atomic configurationthat is infinitely periodic in the xy plane and isolated inthe z direction. The electronic motion is described usingthe Bloch orbital in the slab approximation, u slab n k ( r , t ),where k is the 2D crystal wave vector. This satisfies thefollowing TDKS equation: i (cid:126) ∂u slab n k ( r , t ) ∂t = (cid:26) m (cid:16) − i (cid:126) ∇ + (cid:126) k + ec A ( r , t ) (cid:17) − eφ ( r , t ) + δV ion ( r ) + V xc ( r , t ) (cid:27) u slab n k ( r , t ) . (5)The vector potential A ( r , t ) and the scalar potential φ ( r , t ) are 2D periodic and satisfy the Maxwell equationsin the Coulomb gauge: (cid:18) c ∂ ∂t − ∇ (cid:19) A ( r , t ) + 1 c ∂∂t ∇ φ ( r , t ) = 4 πc j ( r , t ) , (6) ∇ φ ( r , t ) = − πρ ( r , t ) . (7)The vector potential A ( r , t ) is smoothly connected to theasymptotic field of Eq. (1) at planes that are sufficientlyseparated from the medium. The scalar potential is cho-sen to be periodic in the z direction as well as in the x and y directions. As mentioned above, all the dynamicvariables, Bloch orbitals, and scalar and vector poten-tials are calculated using a common spacial grid withoutcoarse graining. In this sense, both the vector and scalarpotentials are treated microscopically in this approach.The charge density ρ ( r , t ) and the electric current den-sity j ( r , t ) are derived from the Bloch orbitals, as follows: ρ ( r , t ) = ρ ion ( r ) − e occ (cid:88) n k | u slab n k ( r , t ) | , (8) j ( r , t ) = − em Re occ (cid:88) n k u slab ∗ n k ( r , t ) (cid:16) − i (cid:126) ∇ + (cid:126) k + ec A ( r , t ) (cid:17) u slab n k ( r , t ) , (9)where ρ ion ( r ) is the charge density of the ion cores. Inthis method, we ignore the current from the nonlocal partof the pseudopotential .Because the Bloch orbitals { u slab n k ( r , t ) } are defined inthe calculation box that includes entire film as well asthe vacuum region, the size of the computation rapidlygrows as the thickness of the film increases. In terms ofthe computational cost, applications of the single-scaleMaxwell-TDDFT method is limited to thin films of thick-ness less than several tens of nanometers. D. Multiscale Maxwell-TDDFT method
In the multiscale Maxwell-TDDFT method , macro-scopic light propagation and microscopic electronic mo-tion are coupled using a coarse-graining approximation.The light propagation is described using the followingmacroscopic wave equation: (cid:18) c ∂ ∂t − ∂ ∂Z (cid:19) A Z ( t ) = 4 πc J Z ( t ) , (10)where Z is the macroscopic coordinate. This wave equa-tion is solved using a one-dimensional grid. At each gridpoint of Z , a bulk system with 3D periodicity is consid-ered. The electronic motion at each grid point of Z isdescribed using 3D-periodic Bloch orbitals, u bulk n k ,Z ( r , t ),which satisfy the TDKS equation: i (cid:126) ∂u bulk n k ,Z ( r , t ) ∂t = (cid:26) m (cid:16) − i (cid:126) ∇ + (cid:126) k + ec ˆ x A Z ( t ) (cid:17) − eφ Z ( r , t ) + δV ion ( r ) + V xc , Z ( r , t ) (cid:27) u bulk n k ,Z ( r , t ) . (11)From the Bloch orbitals, the averaged electric currentdensity J Z ( t ) is calculated in the same manner as Eq. (3),but replacing A (i) ( t ) with A Z ( t ).At the beginning of the calculation, the Bloch or-bitals at each grid point Z , u bulk n k ,Z ( r , t ), are set to theground state. The vector potential A Z ( t ) is set to in-clude only the incident pulse in the vacuum region. Bysolving Eqs. (10) and (11) simultaneously, we can evolveboth A Z ( t ) and u bulk n k ,Z ( r , t ) simultaneously. We notethat microscopic electronic systems at different Z po-sitions interact only through the vector potential A Z ( t ).When the light pulse is sufficiently weak and thus theperturbation approximation is applicable, the multiscale Maxwell-TDDFT method results in the ordinary macro-scopic electromagnetism with the constitutive relationgiven by the dielectric function in the linear responseTDDFT . The multiscale Maxwell-TDDFT method hasbeen successfully applied to investigate a number of ex-tremely nonlinear and ultrafast phenomena including at-tosecond science , saturable absorption , coherentphonon generation and detection , and initial stageof nonthermal laser processing .In the multiscale Maxwell-TDDFT method, the ef-fect of the light propagation can be taken into accountwhile the effect of the surface electronic structure can-not. In compensation for the surface effect, the multi-scale Maxwell-TDDFT method can be applied to thickfilms for which the single-scale Maxwell-TDDFT methodis not capable. As we will show later, results of the multi-scale Maxwell-TDDFT method coincide reasonably withthose of the single-scale Maxwell-TDDFT method whenthe surface effect is not significant. E. 2D approximation
For sufficiently thin films, we may assume that themacroscopic electric field is spatially uniform inside thethin film. We call this the 2D approximation. Thisapproximation is useful to identify and distinguish theconductive effect from the propagation effect, and tounderstand the thickness dependence of the HHG. The2D approximation can be derived from the multiscaleMaxwell-TDDFT method, as described below. Alter-natively, it can also be derived from the single-scaleMaxwell-TDDFT method, as explained in Ref 47.We consider a film that is sufficiently thick to justifyneglecting the effect of the surface electronic structureand sufficiently thin to regard the macroscopic electricfield as spatially uniform in the film. In the next section,we see that there are thickness regions in which bothassumptions are fulfilled simultaneously. Under these as-sumptions, the current density in Eq. (10) may be treatedas J Z ( t ) (cid:39) δ ( Z ) J ( t ) d, (12)where d is the thickness of the film and J ( t ) is the currentdensity averaged over the unit cell. By inserting thiscurrent density into Eq. (10), we have (cid:18) c ∂ ∂t − ∂ ∂Z (cid:19) A Z ( t ) = 4 πc δ ( Z ) J ( t ) d. (13)This equation should be solved with the asymptoticform of the vector potential of Eq. (1). In practice, thevector potential of Eq. (1) is the solution of Eq. (13),except at Z = 0. The reflected and transmitted fieldsare determined by the connection conditions at Z = 0.From the continuity of the vector potential at Z = 0, wehave A (i) ( t ) + A (r) ( t ) = A (t) ( t ) . (14)By integrating Eq. (13) once over Z , we obtain dA (t) dt = dA (i) dt + 2 πdJ [ A (t) ]( t ) (15)at Z = 0, where we denote the current density as J [ A (t) ]( t ) to indicate that this is the current densitycaused by the vector potential A Z =0 ( t ) = A (t) ( t ). Thisis the basic equation of the 2D approximation.We note that this 2D approximation is equivalent to aspecific case of the multiscale Maxwell-TDDFT method.If we take only a single grid point in the macroscopiccoordinate Z for the medium and choose the grid spacingto be equal to d , the calculation using the multiscaleMaxwell-TDDFT method coincides with the above 2Dapproximation.From the 2D approximation, we obtain several use-ful understandings regarding the HHG from very thinfilms. The continuity equation (14) indicates that theHHG spectra of the reflected and transmitted pulses areequal. Equation (15) indicates that the current densitythat creates HHG is produced by the transmitted pulse A (t) ( t ), not the incident pulse of A (i) ( t ), and the trans-mitted pulse itself is modulated by the current density.Even in very thin films such as monatomic layered films,the transmitted (and simultaneously, the reflected) fieldsare modulated by the current that flows in the film. Thismodulation significantly affects the HHG spectrum, asshown in the next section. We call this effect the conduc-tive effect, as it is related to the electric current flowingin the thin film. F. High harmonic generation spectrum
In the single-scale and multiscale Maxwell-TDDFTmethods, we directly obtain the vector potentials of thereflected and the transmitted pulses. In the 2D approx-imation, we obtain the vector potentials of the trans-mitted pulse by Eq. (15). The vector potential of thereflected pulse is obtained by using the relation Eq. (14).From the vector potentials, we obtain the respective elec-tric fields by taking the time derivative.For the single unit-cell method, we will use the follow-ing expression for the reflected and transmitted electricfields. E (r) ( t ) = E (t) ( t ) − E (i) ( t ) = − πdc J [ A (i) ]( t ) . (16)This relation is derived if we replace A (t) with A (i) in J [ A (t) ]( t ) of Eq. (15).We evaluate the HHG spectra via the square ofthe Fourier-transformed electric field | E (r , t) ( ω ) | , where E (r , t) ( ω ) is defined by, E (r , t) ( ω ) = (cid:90) t + Tt dt e iωt E (r , t) ( t ) f (cid:18) t − t T (cid:19) . (17) Here, f ( x ) ≡ − x + 2 x is a smoothing function and t is the initial time. For the reflection pulse, t is set tothe initial time of the incident pulse. For the transmittedpulse, the time delay by transmission is added to t forthick films. When we discuss the thickness dependenceof HHG, we introduce the strength of the n th-order har-monics, as follows: I (r , t) n = (cid:90) ( n + ) ω ( n − ) ω dω | E (r , t) ( ω ) | , (18)where ω is the fundamental frequency. G. Numerical detail
For numerical calculations, we utilize the open-sourcesoftware SALMON (Scalable Ab initio Light-Matter sim-ulator for Optics and Nanoscience) developed by ourgroup. In this code, electronic orbitals as well as electro-magnetic fields are expressed using a uniform 3D spatialgrid. The time evolution of the electron orbitals is carriedout using the Taylor expansion method .In the single unit-cell method, we use a cubic unit cellthat contains eight Si atoms with the side length of a =0 .
543 nm. To express the Bloch orbitals, a spatial grid of16 points is used. The 3D Brillouin zone is sampled bya 24 k -point grid. The time step is set to 2 . × − fs.The same cubic unit cell is also used in the 2D ap-proximation and multiscale Maxwell-TDDFT method.In the multiscale Maxwell-TDDFT method, a uniform1D grid is also introduced to describe the wave equationof Eq. (10), with the grid spacing less than or equal to6 .
25 nm.In the single-scale Maxwell-TDDFT method, the 2D-periodic box of the size of a × a × ( d +4 a ) is used where thesize along the z axis contains the vacuum region of 4 a inthe slab approximation. The atomic positions of Si atomsare set at those positions in the bulk crystalline system,and the dangling bonds at the surfaces are terminated byhydrogen atoms. A uniform 3D spatial grid is used withthe same grid spacing as that in the cubic unit cell. The2D Brillouin zone of the single-scale Maxwell-TDDFTmethod is sampled by a 24 × k -point grid. III. RESULTS AND DISCUSSIONA. HHG in reflected and transmitted pulses
We employ the following time profile for the incidentpulse: A (i) ( t ) = − cE ω sin ω t cos (cid:18) πtT (cid:19) , ( − T < t < T , (19)where E is the maximum amplitude of the electric field, ω is the average frequency, and T is the pulse duration. z (μm) incident pulsereflected pulse transmitted pulse (b)(c) incident spectra harmonic order → (a) x zy Si film (d = 10a)incident pulse −0.3−0.2−0.1 0 0.1 0.2 0.3 A / c ( a . u . ) −0.2−0.1 0 0.1 0.2 0.3 −8 −6 −4 −2 0 2 4 6 8 A / c ( a . u . )
1 3 5 7 9 11 13 15 17 1 3 5 7 9 11 13 15 17
RHHG THHG
1 3 5 7 9 11 13 15 17 harmonic order harmonic order
FIG. 1. Typical calculation using the single-scale Maxwell-TDDFT method for a Si thin film of thickness d = 10 a (5.43nm). (a) Electron density in the ground state in a planecontaining Si atoms. (b) Initial pulse set to the left of the thinfilm located at z = 0. The peak intensity of the pulse is set at I = 5 × W/cm at t = 0. (c) Reflected and transmittedpulses at t = 31 .
25 fs. The insets show the correspondingHHG spectra.
In the following calculations, we set (cid:126) ω = 1.5 eV and T = 30 fs.For a typical case, Figure 1 shows snapshots of theincident pulse [Eq. (19)] and the transmitted and thereflected pulses in the calculation using the single-scaleMaxwell-TDDFT method for a film of thickness d = 10 a (5.43 nm). The maximum amplitude E of the incidentpulse is set to provide the peak intensity of I = 5 × W/cm . In Fig. 1(a), the electronic density at theground state in a plane that includes Si atoms is shown.In Fig. 1(b), the incident pulse is prepared to the left ofthe thin film located at z = 0. Figure 1(c) shows thepulses at a sufficient time after interaction. The reflected(transmitted) pulse is seen in the left (right) to the film.The insets show spectra for the respective pulses. It canbe clearly seen that the reflected and the transmittedpulses include HHG signals.We first discuss thin films with thicknesses less thana few tens of nanometers. For such films, it is possible to achieve calculations using the single-scale Maxwell-TDDFT method. Figure 2 shows the HHG spectra in-cluded in the reflected (RHHG) and transmitted (THHG)pulses. The left panels [Fig. 2(a)–(d)] show the spectrafor the incident pulse with a peak intensity of I = 1 × W/cm , and the right panels [Fig. 2(e)–(h)] show thosefor a peak intensity of 5 × W/cm . Panels (a) and(e) show results for the film of thickness d = a (0.543nm), (b) and (f) for 10 a (5.43 nm), (c) and (g) for 30 a (16.29 nm), and (d) and (h) for 50 a (27.15 nm). Thered solid line (blue dashed line) shows RHHG (THHG)using the single-scale Maxwell-TDDFT method, and theblack dotted line shows HHG using the 2D approximation[Eq. (15)]. In the 2D approximation, the HHG spectrain the reflected and transmitted pulses coincide, exceptfor the component of the fundamental frequency ω . Theorange dash-dotted line shows the HHG spectrum usingthe single unit-cell method.First, we focus on the extremely thin case of d = a (0.543 nm) [Fig. 2(a), (e)]. This is a film composed of twoatomic layers. Although fabrication of such thin materialcomposed of Si is not easy, HHG of such thin materialsis highly concerned, as measurements and calculationshave been reported for a number of 2D materials suchas graphene and transition metal dichalcogenides .Here the RHHG and the THHG spectra using the single-scale Maxwell-TDDFT method coincide accurately, asexpected from the argument in the 2D approximation[Eq. (14)]. However, the HHG spectrum by the 2D ap-proximation (black dotted line) does not coincide accu-rately with that using the single-scale Maxwell-TDDFTmethod. This indicates the significance of the effect ofthe surface electronic structure that is not included inthe 2D approximation. We note that the validity of therelation of Eq. (14) does not suffer by the surface elec-tronic structure, as it can be derived directly from thesingle-scale Maxwell-TDDFT method . The result ofthe single unit-cell method (orange dash-dotted line) co-incides with the 2D approximation, indicating that theconductive effect in the film is not important in such ex-tremely thin film.The above observations hold also at the film of thick-ness d = 10 a (5.43 nm) [Fig. 2(b), (f)], except that thediscrepancy between spectra using the 2D approximationand the single unit-cell method becomes apparent. It im-plies that the effect of the conductivity in the thin film,use of A (t) , not the incident pulse A (i) , to evaluate thecurrent in the right hand side of Eq. (15), increases as thefilm thickness increases, because the effect is proportionalto d .In the case of d = 30 a (16.29 nm) [Fig. 2(c), (g)], theresults using the single-scale Maxwell-TDDFT methodand those using the 2D approximation agree reasonablywith each other. This indicates that neither the surfaceelectronic structure nor the propagation effect is signifi-cant at this thickness. Only the conductive effect is im-portant. By closely observing the difference between theRHHG and THHG spectra, the latter is slightly greater -18 -15 -12 -9 -6 -3 | E ( ω ) | ( a . u . ) (a) Single-scale Maxwell-TDDFT (RHHG)Single-scale Maxwell-TDDFT (THHG)2D approximationSingle unit-cell -18 -15 -12 -9 -6 -3 | E ( ω ) | ( a . u . ) (b)10 -18 -15 -12 -9 -6 -3 | E ( ω ) | ( a . u . ) (c)10 -18 -15 -12 -9 -6 -3
1 3 5 7 9 11 13 15 17 | E ( ω ) | ( a . u . ) Harmonic order(d) 10 -15 -12 -9 -6 -3 (e)10 -15 -12 -9 -6 -3 (f)10 -15 -12 -9 -6 -3 (g)10 -15 -12 -9 -6 -3
1 3 5 7 9 11 13 15 17Harmonic order(h)Peak intensity I = 1×10 W/cm Peak intensity I = 5×10 W/cm d = ad = 10ad = 30ad = 50a FIG. 2. HHG spectra in the reflected pulse (RHHG) and the transmitted pulse (THHG). Left panels [(a)–(d)] show the spectrafor a pulse with the peak intensity of I = 1 × W/cm and right panels [(e)–(h)] for a pulse of I = 5 × W/cm . Toppanels [(a), (e)] show the spectra for the film thickness of d = a (0.543 nm), second panels [(b), (f)] for d = 10 a (5.43 nm),third panels [(c), (g)] for d = 30 a (16.29 nm), and bottom panels [(d), (h)] for d = 50 a (27.15 nm). The red solid (blue dashed)line for the RHHG (THHG) using the single-scale Maxwell-TDDFT method, blue dotted line using the 2D approximation, andorange dash-dotted line using the single unit-cell method. than the former. This indicates that the propagation ef-fect starts to appear.In the case of d = 50 a (27.15 nm) [Fig. 2(d), (h)], it isclear that, in the single-scale Maxwell-TDDFT method,the THHG is greater than the RHHG due to the propa-gation effect. Now the calculation using the single unit-cell method greatly differs from the others, despite the2D approximation still reasonably reproduces the resultsusing the single-scale Maxwell-TDDFT method. Look-ing into detail, the spectrum of the 2D approximation iscloser to the THHG in the single-scale Maxwell-TDDFTmethod than the RHHG. The high-order part of the spec-trum shown in Fig. 2(d) [and partially Fig. 2(c)] looksnoisy, presumably because the driving field is too weak( I = 1 × W/cm ) to produce high-order signals atthis thickness.Overall, the 2D approximation is a reasonable approx-imation to the single-scale Maxwell-TDDFT method inthis thickness region. This indicates that the surface elec-tronic structure and the propagation effects are not sig-nificant. The large difference between the results using the 2D approximation and those using the single unit-cellmethod indicates that the conductive effect is significant,even in these very thin films.By comparing the different thicknesses, especiallyFig. 2(e) and (f)–(h), we realize that the HHG spectrumis unclear in films of thickness d = a (0.543 nm) and be-comes clearer as the thickness increases. This trend iscommon in both spectra using the single-scale Maxwell-TDDFT method and the 2D approximation. It impliesthat the conductive effect, the current in the thin filmthat appears in the right hand side of Eq. (15), con-tributes to produce a cleaner HHG signal. This is con-sistent with the observation in Ref 35 where both theconductive and the propagation effects are included. B. Thickness dependence
Figure 3 shows the thickness dependence of the inten-sity of HHG in the reflected pulse using Eq. (18). Thepeak intensity of the incident pulses is set at I = 1 × -7 -6 -5 I ( a . u . ) (a) Single-scale Maxwell-TDDFT2D approximation -10 -9 -8 I ( a . u . ) (b)10 -13 -12 -11 I ( a . u . ) (c)10 -16 -15 -14 -13 -12
0 5 10 15 20 25 30 35 40 I ( a . u . ) Thickness (nm)(d) 10 -6 -5 -4 I ( a . u . ) (e) Single-scale Maxwell-TDDFT2D approximation -7 -6 -5 I ( a . u . ) (f)10 -9 -8 -7 I ( a . u . ) (g)10 -10 -9
0 5 10 15 20 25 30 35 40 I ( a . u . ) Thickness (nm)(h)Peak intensity I = 1×10 W/cm Peak intensity I = 5×10 W/cm FIG. 3. Thickness dependence of HHG intensities is shown from 3rd [(a), (e)], 5th [(b), (f)], 7th [(c), (g)], and 9th [(d), (h)]orders included in the reflected pulse. Left panels [(a)–(d)] for the incident pulse of the peak intensity at I = 1 × W/cm ,and right panels [(e)–(h)] at I = 5 × W/cm . Results using single-scale Maxwell-TDDFT method are shown by red solidlines, and results using 2D approximation are shown by blue-dashed lines. W/cm in Fig. 3(a)–(d) and at 5 × W/cm inFig. 3(e)–(h). Intensities of harmonic order of 3rd [(a),(e)], 5th [(b), (f)], 7th [(c), (g)], and 9th [(d), (h)] areshown. The red solid line (blue dashed line) correspondsto the single-scale Maxwell-TDDFT method (2D approx-imation).As seen from Fig. 3, the thickness dependence of HHGintensities using the 2D approximation (blue dashed line)shows qualitative agreement with those using the single-scale Maxwell-TDDFT method (red solid line). In bothcalculations, the HHG intensity is maximum aroundthicknesses of d = 2–15 nm, irrespective of the orderof the HHG. The appearance of the maximum at thesethicknesses can be understood as a consequence of theconductive effect, as described below, extending the for-malism of the 2D approximation.By taking the Fourier transform of Eq. (15), we obtain − iωA (t) ( ω ) = − iωA (i) ( ω ) + 2 πd J [ A (t) ]( ω ) . (20)We decompose the current density into linear andnonlinear components as J [ A (t) ]( ω ) = J L [ A (t) ]( ω ) + J NL [ A (t) ]( ω ) and use the constitutive relation for the lin- ear part, J L [ A (t) ]( ω ) = σ ( ω ) E (t) ( ω ) = iωc σ ( ω ) A (t) ( ω ) , (21)where σ ( ω ) is the conductivity of the bulk medium. Com-bining the above equations, we obtain (cid:18) πdc σ ( ω ) (cid:19) A (t) ( ω ) = A (i) ( ω ) + 2 πidω J NL [ A (t) ]( ω ) . (22)We also decompose the transmitted vector potential intolinear and nonlinear components , A (t) ( ω ) = A (t)L ( ω ) + A (t)NL ( ω ), where the linear vector potential is given by A (t)L ( ω ) = (cid:18) πdc σ ( ω ) (cid:19) − A (i) ( ω ) . (23)The nonlinear term A (t)NL satisfies the following relation: A (t)NL ( ω ) = 2 πidω (cid:18) πdc σ ( ω ) (cid:19) − J NL [ A (t)L + A (t)NL ]( ω ) . (24)We note that the nonlinear component of the reflectedpulse is also given by Eq. (24). So far, no approximationshave been made besides the 2D approximation. Here, weassume that the nonlinear component A (t)NL is rather smalland can be ignored to evaluate the nonlinear current den-sity in the right hand side of Eq. (24). Then, we have A (t)NL ( ω ) (cid:39) πidω (cid:18) πdc σ ( ω ) (cid:19) − J NL [ A (t)L ]( ω ) . (25)To investigate how the HHG intensity included in theabove A (t)NL ( ω ) depends on d , we introduce two furtherassumptions. First, we assume that the incident pulse A (i) ( t ) has a well-defined frequency, which we denote as ω . Second, the n th-order nonlinear current density isproportional to A n , where A is the amplitude of thevector potential. Then, from Eq. (23), the d depen-dence of the n th-order term in J NL [ A (t)L ]( ω ) is given by(1 + 2 πσ ( ω ) d/c ) − n . Then, from Eq. (25), the n th-orderterm of A (t)NL ( ω = nω ) is proportional to the followingfactor: A (t)NL ( nω ) ∝ d (cid:16) πσ ( nω ) c d (cid:17) (cid:16) πσ ( ω ) c d (cid:17) n . (26)From this d dependence, we expect that the intensity ofthe HHG of each order shows a maximum as a functionof d , and the peak position is determined by the bulkconductivity.The conductivity at the frequency of the HHG, σ ( nω ),depends on the order n and becomes small at high orders, n >>
1. For simplicity, we consider the d -dependencecaused by the conductivity at the fundamental frequency, σ ( ω ). We expect the peak of HHG appears at the thick-ness, d ∼ c π √ n − | σ ( ω ) | , (27)where we used σ ( ω ) as a pure imaginary number at thefrequency (cid:126) ω = 1 . c/ (2 π | σ ( ω ) | ) (cid:39) . ω . This explains the appearance of the peak around d =2–15 nm. C. Propagation effect
As the film thickness increases, the computational costof the single-scale Maxwell-TDDFT method rapidly in-creases. Therefore, the use of the multiscale Maxwell-TDDFT method is appropriate and necessary. To con-firm the accuracy of the multiscale method, we comparethe two methods in Fig. 4 for RHHG and THHG emittedfrom a film of thickness d = 50 a (27.19 nm) and the inci-dent pulse with maximum intensity I = 5 × W/cm ,which is the same conditions as those in Fig. 2(h). Be-cause the spectra using two methods coincide accurately -15 -12 -9 -6 -3 | E ( ω ) | ( a . u . ) (a) RHHG Single-scale Maxwell-TDDFTMultiscale Maxwell-TDDFT10 -15 -12 -9 -6 -3
1 3 5 7 9 11 13 15 17 | E ( ω ) | ( a . u . ) Harmonic order(b) THHG Single-scale Maxwell-TDDFTMultiscale Maxwell-TDDFT
FIG. 4. HHG spectra of d = 50 a (27.19 nm) film includedin the reflected (a) and the transmitted (b) pulses are shown.Red solid lines show the results using the single-scale Maxwell-TDDFT method, whereas the blue dashed lines show thoseusing the multiscale Maxwell-TDDFT method. with each other, we may conclude that the multiscaleMaxwell-TDDFT method, which ignores the effect of thesurface electronic structure, is reliable for films of thisthickness and greater.We note that, as mentioned in Sec. II.C, the multi-scale Maxwell-TDDFT method is identical to the 2D ap-proximation for very thin films where the macroscopicelectric field may be regarded as uniform inside the thinfilm. As discussed in Fig. 2, the 2D approximation wasa good approximation for films of thickness less than d = 50 a (27.19 nm). In summary, we expect the mul-tiscale Maxwell-TDDFT method to be reliable for thinfilms of any thickness as long as the surface electronicstructure is not important.Figure 5(a)–(d) shows the thickness dependence of thethird- to ninth-order harmonics included in the reflectedand transmitted pulses for the incident pulse with a peakintensity of I = 5 × W/cm . Results using the single-scale Maxwell-TDDFT method are shown for thickness d ≤ a (27.19 nm). The RHHG shown here are the sameas those in Fig. 3(e)–(h). The results using the multiscaleMaxwell-TDDFT method are shown for thickness d ≥ a (16.29 nm).Figure 5(e) shows the square of the maximum electricfield amplitude in time, E max = max t | E ( t ) | , at the frontsurface and the back surface of the film. When using thesingle-scale Maxwell-TDDFT method, an average overthe surface is taken. Figure 5(f) shows | t | and | r | , where t and r are the amplitude of transmission andreflection coefficients in ordinary electromagnetism: r = r (1 − e iφ )1 − r e iφ , t = (1 − r ) e iφ − r e iφ , (28)where r = (1 − n ) / (1 + n ) and φ = 2 πnd/λ . The index0 -7 -6 -5 -4 I ( a . u . ) (a)10 -9 -7 -5 I ( a . u . ) (b)10 -11 -9 -7 I ( a . u . ) (c)10 -13 -11 -9 I ( a . u . ) (d)10 -2 -1 E m a x2 ( a . u . ) (e) front surface (single-scale )front surface (multiscale)back surface (single-scale)back surface (multiscale) 0 0.2 0.4 0.6 0.8 1 0 50 100 150 200 R a t e Thickness (nm)(f) T=|t| |1+r| RHHG (single-scale)RHHG (multiscale)THHG (single-scale)THHG (multiscale)
FIG. 5. (a)–(d) Thickness dependence of the intensities ofthe reflected (solid lines with crosses) and the transmitted(dashed lines with squares) HHG, from 3rd (a) to 9th (d)order. Red lines and symbols are calculated using the single-scale Maxwell-TDDFT method, and blue lines and symbolsusing the multiscale Maxwell-TDDFT method. (e) Maximumintensities of the electric fields averaged over the 2D area ofthe front (back) surface are shown by red lines with sym-bols (blue lines with symbols). (f) Intensities of electric fieldsin ordinary electromagnetism at the front (back) surface areshown by dashed (solid) lines. of refraction n of Si and the wavelength λ is evaluated atthe frequency ω . T = | t | is equal to the transmittanceand is equal to the square of the electric field at the backsurface, while | r | is equal to the square of the electricfield at the front surface.Because the frequency ω is below the direct bandgap of Si, | t | and | r | become unity when 2 πnd is equal tothe wavelength and its integer multiples. As seen fromFig. 5(e), however, the magnitude of the electric fieldat the back surface decreases monotonically as the filmthickness increases in the Maxwell-TDDFT calculations.The magnitude of the electric field at the front surfaceshows more oscillatory behavior, with a clear minimumat πd = λ . However, the magnitude is smaller at πd = 2 λ than at the zero-thickness limit. These decreases origi-nate the nonlinear propagation effect included in boththe single-scale and multiscale Maxwell-TDDFT meth-ods. In any case, the magnitude of the electric field atneither the front nor the back surface exceeds the mag-nitude of the electric field at the zero-thickness limit.The intensities of RHHG and THHG coincide witheach other when the film is sufficiently thin. As discussedpreviously, they are well understood by Eq. (14) and startto deviate at a distance of 10–20 nm. As seen from Fig. 5,except the d <
20 nm region, the thickness dependenceof RHHG (THHG) is correlated faithfully with the max-imum electric field at the front (back) surface. This indi-cates that the magnitude of the RHHG and the THHG isdetermined simply by the intensity of its driving electricfield at the respective surfaces. Because the magnitudeof the electric field at the front or back surface does notexceed those at the zero-thickness limit, the maximumof the intensity of the HHG signals also appears around d ∼
10 nm as seen in Fig. 3.
IV. CONCLUSION
We have developed a few theoretical and computa-tional methods based on first-principles TDDFT and in-vestigated HHG from thin films of crystalline solids withthicknesses from a few atomic layers to a few hundredsof nanometers. Using these methods, it is possible toinvestigate effects of (1) surface electronic structure, (2)conductivity of the film, and (3) light propagation of fun-damental and high-harmonic fields, as well as (4) the ef-fect of electronic motion in the bulk energy band.Among methods used herein, the most sophisticatedmethod is the single-scale Maxwell-TDDFT method,which considers all four effects mentioned above. It canbe applied to thin films with thicknesses less than a fewtens of nanometers. For thicker films, the multiscaleMaxwell-TDDFT method, in which a coarse-graining ap-proximation is introduced, is applicable. It can treatthe effects except for the surface electronic structure.To achieve insight into the thickness dependence of theHHG, the 2D approximation is useful. We applied thesemethods to thin films of crystalline silicon as a prototypematerial, and obtained the following conclusions.For thin films with thicknesses less than a few tensof nanometers, HHG spectra of reflected and transmit-ted pulses are almost identical. This could be under-stood from the 2D approximation. In extremely thinfilms, the HHG intensity increases as the thickness of1the film increases. However, the intensity soon saturatesand reaches its maximum at a thickness of around 2–15nm. The saturation of the intensity originates from thecurrent that flows in the thin film and can be describedusing the bulk linear conductivity of the medium. Thisconductive effect also creates a clean HHG spectrum.As the thickness increases beyond a few tens ofnanometers, it is found that the intensity of the reflected(transmitted) HHG strongly correlates with the intensityof the electric field at the front (back) surface of the thinfilm. By reflecting the interference effect in pulse prop-agation, we find an enhancement in the reflected HHGwhen the film thickness is equal to integer or half-integermultiples of the fundamental wavelength in the medium, λ/n , where n is the index of refraction. However, thetransmitted HHG monotonically decreases as the thick- ness increases, owing to the nonlinear propagation of thefundamental wave. ACKNOWLEDGMENTS
This research was supported by JST-CREST undergrant number JP-MJCR16N5, and by MEXT Quan-tum Leap Flagship Program (MEXT Q-LEAP) GrantNumber JPMXS0118068681, and by JSPS KAKENHIGrant Number 20H2649. Calculations are carried outat Oakforest-PACS at JCAHPC under the supportthrough the HPCI System Research Project (Project ID:hp20034), and Multidisciplinary Cooperative ResearchProgram in CCS, University of Tsukuba. P. B. Corkum, Phys. Rev. Lett. , 1994 (1993). M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, andP. B. Corkum, Phys. Rev. A , 2117 (1994). T. Brabec and F. Krausz, Rev. Mod. Phys. , 545 (2000). F. Krausz and M. Ivanov, Rev. Mod. Phys. , 163 (2009). S. Ghimire and D. A. Reis, Nat. Phys. , 10 (2019). J. Li, J. Lu, A. Chew, S. Han, J. Li, Y. Wu, H. Wang,S. Ghimire, and Z. Chang, Nat. Commun. , 2748 (2020). S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F.DiMauro, and D. A. Reis, Nat. Phys. , 138 (2011). O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek,C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W.Koch, and R. Huber, Nat. Photonics , 119 (2014). G. Vampa, T. J. Hammond, N. Thir´e, B. E. Schmidt,F. L´egar´e, C. R. McDonald, T. Brabec, and P. B. Corkum,Nature , 462 (2015). M. Hohenleutner, F. Langer, O. Schubert, M. Knorr,U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature , 572 (2015). T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T.Hassan, and E. Goulielmakis, Nature , 498 (2015). F. Langer, M. Hohenleutner, U. Huttner, S. W. Koch,M. Kira, and R. Huber, Nat. Photonics , 227 (2017). Y. S. You, D. A. Reis, and S. Ghimire, Nat. Phys. , 345(2017). H. Shirai, F. Kumaki, Y. Nomura, and T. Fuji, Opt. Lett. , 2094 (2018). G. Vampa, T. J. Hammond, M. Taucer, X. Ding, X. Ropag-nol, T. Ozaki, S. Delprat, M. Chaker, N. Thir´e, B. E.Schmidt, F. L´egar´e, D. D. Klug, A. Y. Naumov, D. M.Villeneuve, A. Staudte, and P. B. Corkum, Nat. Photon-ics , 465 (2018). G. Orenstein, A. Julie Uzan, S. Gadasi, T. Arusi-Parpar,M. Kr¨uger, R. Cireasa, B. D. Bruner, and N. Dudovich,Opt. Express , 37835 (2019). N. Yoshikawa, T. Tamaya, and K. Tanaka, Science (80-.). , 736 (2017). N. Yoshikawa, K. Nagai, K. Uchida, Y. Takaguchi,S. Sasaki, Y. Miyata, and K. Tanaka, Nat. Commun. ,3709 (2019). H. Liu, C. Guo, G. Vampa, J. L. Zhang, T. Sarmiento,M. Xiao, P. H. Bucksbaum, J. Vuˇckovi´c, S. Fan, and D. A. Reis, Nat. Phys. , 1006 (2018). C. Yu, S. Jiang, and R. Lu, Adv. Phys. X , 1562982(2019). K. K. Hansen, T. Deffge, and D. Bauer, Phys. Rev. A ,053418 (2017). K. K. Hansen, D. Bauer, and L. B. Madsen, Phys. Rev.A , 043424 (2018). J.-Z. Jin, H. Liang, X.-R. Xiao, M.-X. Wang, S.-G. Chen,X.-Y. Wu, Q. Gong, and L.-Y. Peng, Phys. Rev. A ,013412 (2019). T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto,M. Kuwata-Gonokami, and K. L. Ishikawa, Phys. Rev.A , 043416 (2017). M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B.Gaarde, Phys. Rev. A , 043839 (2015). T. Apostolova and B. Obreshkov, Diam. Relat. Mater. ,165 (2018). S. Ghimire, A. D. DiChiara, E. Sistrunk, G. Ndabashimiye,U. B. Szafruga, A. Mohammad, P. Agostini, L. F. Di-Mauro, and D. A. Reis, Phys. Rev. A , 043836 (2012). G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B.Corkum, and T. Brabec, Phys. Rev. Lett. , 073901(2014). G. Vampa, C. R. McDonald, G. Orlando, P. B. Corkum,and T. Brabec, Phys. Rev. B , 064302 (2015). F. H. M. Faisal and J. Z. Kami´nski, Phys. Rev. A , 748(1997). F. H. M. Faisal, J. Z. Kami´nski, and E. Saczuk, Phys.Rev. A , 023412 (2005). T. Otobe, J. Appl. Phys. , 093112 (2012). T. Otobe, Phys. Rev. B , 235152 (2016). N. Tancogne-Dejean, O. D. M¨ucke, F. X. K¨artner,and A. Rubio, Phys. Rev. Lett. , 087403 (2017),arXiv:1609.09298. I. Floss, C. Lemell, G. Wachter, V. Smejkal, S. A. Sato,X.-M. Tong, K. Yabana, and J. Burgd¨orfer, Phys. Rev. A , 011401 (2018). A. T. Georges and N. E. Karatzas, Appl. Phys. B , 479(2005). G. Le Breton, A. Rubio, and N. Tancogne-Dejean, Phys.Rev. B , 165308 (2018). O. Neufeld, D. Podolsky, and O. Cohen, Nat. Commun. , 405 (2019). T. Popmintchev, M.-C. Chen, A. Bahabad, M. Ger-rity, P. Sidorenko, O. Cohen, I. P. Christov, M. M.Murnane, and H. C. Kapteyn, Proceedings of theNational Academy of Sciences G. Vampa, Y. S. You, H. Liu, S. Ghimire, and D. A. Reis,Opt. Express , 12210 (2018). P. Xia, C. Kim, F. Lu, T. Kanai, H. Akiyama, J. Itatani,and N. Ishii, Opt. Express , 29393 (2018). P. A. Zhokhov and A. M. Zheltikov, Phys. Rev. A ,033415 (2017). G. Orlando, T.-S. Ho, and S.-I. Chu, J. Opt. Soc. Am. B , 1873 (2019). I. Kilen, M. Kolesik, J. Hader, J. V. Moloney, U. Huttner,M. K. Hagen, and S. W. Koch, Phys. Rev. Lett. ,083901 (2020). E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984). C. A. Ullrich,
Time-Dependent Density-Functional Theory:Concepts and Applications (Oxford Graduate Texts) (Ox-ford Univ Pr (Txt), 2012). S. Yamada, M. Noda, K. Nobusada, and K. Yabana, Phys.Rev. B , 245147 (2018). K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe,and G. F. Bertsch, Phys. Rev. B , 045134 (2012),arXiv:arXiv:1112.2291v1. G. F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, Phys.Rev. B , 7998 (2000), arXiv:0005512v1 [arXiv:cond-mat]. T. Otobe, M. Yamagiwa, J.-I. Iwata, K. Yabana, T. Nakat-sukasa, and G. F. Bertsch, Phys. Rev. B , 165104(2008). M. Uemoto, Y. Kuwabara, S. A. Sato, and K. Yabana,The Journal of Chemical Physics , 094101 (2019),https://doi.org/10.1063/1.5068711. A. Yamada and K. Yabana, The European Physical Jour-nal D , 87 (2019). G. Wachter, C. Lemell, J. Burgd¨orfer, S. A. Sato, X.-M. Tong, and K. Yabana, Phys. Rev. Lett. , 087401(2014). M. Schultze, K. Ramasesha, C. Pemmaraju, S. Sato,D. Whitmore, A. Gandman, J. S. Prell, L. J.Borja, D. Prendergast, K. Yabana, D. M. Neu- mark, and S. R. Leone, Science , 1348 (2014),https://science.sciencemag.org/content/346/6215/1348.full.pdf. S. Yamada and K. Yabana, Phys. Rev. B , 165128(2020). N. Troullier and J. L. Martins, Physical review B , 1993(1991). J. P. Perdew and A. Zunger, Physical Review B , 5048(1981). G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. ,601 (2002). M. Lucchini, S. A. Sato, A. Ludwig, J. Herrmann,M. Volkov, L. Kasmi, Y. Shinohara, K. Yabana,L. Gallmann, and U. Keller, Science , 916 (2016),https://science.sciencemag.org/content/353/6302/916.full.pdf. A. Sommer, E. M. Bothschafter, S. A. Sato, C. Jakubeit,T. Latka, O. Razskazovskaya, H. Fattahi, M. Jobst,W. Schweinberger, V. Shirvanyan, V. S. Yakovlev, R. Kien-berger, K. Yabana, N. Karpowicz, M. Schultze, andF. Krausz, Nature , 86 (2016). M. Uemoto, S. Kurata, N. Kawaguchi, and K. Yabana,Phys. Rev. B , 085433 (2021). A. Yamada and K. Yabana, Phys. Rev. B , 245103(2019). A. Yamada and K. Yabana, Phys. Rev. B , 214313(2020). S. A. Sato, K. Yabana, Y. Shinohara, T. Otobe, K.-M. Lee,and G. F. Bertsch, Phys. Rev. B , 205413 (2015). E. G. Mishchenko, Phys. Rev. Lett. , 246802 (2009),arXiv:0907.3746. K. Yabana, T. Nakatsukasa, J.-I. Iwata, and G. F. Bertsch,Phys. status solidi , 1121 (2006). M. Noda, S. A. Sato, Y. Hirokawa, M. Uemoto,T. Takeuchi, S. Yamada, A. Yamada, Y. Shinohara, M. Ya-maguchi, K. Iida, I. Floss, T. Otobe, K.-M. K.-M. Lee,K. Ishimura, T. Boku, G. F. Bertsch, K. Nobusada, andK. Yabana, Comput. Phys. Commun. , 356 (2019),arXiv:1804.01404. http://salmon-tddft.jp , SALMON official website. K. Yabana and G. F. Bertsch, Phys. Rev. B , 4484(1996). S. A. Sato, H. Hirori, Y. Sanari, Y. Kanemitsu, and A. Ru-bio, Phys. Rev. B103