{'fraction_non_alphanumeric': 0.030487259722843094, 'fraction_numerical': 0.01282968261063925, 'mean_word_length': 4.960831334932054, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Velocity-gauge real-time TDDFT within a numerical atomic orbital basis set

The interaction of laser fields with solid-state systems can be modeled efficiently within the velocity-gauge formalism of real-time time dependent density functional theory (RT-TDDFT). In this article, we discuss the implementation of the velocity-gauge RT-TDDFT equations for electron dynamics within a linear combination of atomic orbitals (LCAO) basis set framework. Numerical results obtained from our LCAO implementation, for the electronic response of periodic systems to both weak and intense laser fields, are compared to those obtained from established real-space grid and Full-Potential Linearized Augumented Planewave approaches. Potential applications of the LCAO based scheme in the context of extreme ultra-violet and soft X-ray spectroscopies involving core-electronic excitations are discussed.


Introduction
Over the last two decades, real-time time dependent density functional theory 1-3 (RT-TDDFT) approaches, wherein the electron density is explicitly propagated in time through numerical integration of the time-dependent Kohn-Sham equations, have gained in prominence as practical first-principles methods for studying electron dynamics in a wide range of quantum systems [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] . By directly simulating density fluctuations in time 2,3 , RT-TDDFT provides a versatile and computationally tractable framework for accessing linear and nonlinear response properties of materials as well as electron dynamics in conditions beyond the perturbative regime, for instance, under the action of intense ultrafast laser pulses 18,19 . While a variety of utilizations of RT-TDDFT have been demonstrated in recent years  , one of the most relevant application domains for RT-TDDFT is the study of laser-matter interactions thanks to the broad appeal of ultrafast laser spectroscopies as experimental tools for investigating and controlling excited states of matter 9,11,13,[18][19][20][21][22][23][24][25] . This is especially the case in a solid-state or condensed-matter context where RT-TDDFT represents perhaps the only computationally feasible first-principles approach for treating the action of intense laser fields 11 . Condensed phase implementations of RT-TDDFT 11,13,26 that rely on the periodic supercell framework treat light-matter interaction through the so-called velocity-gauge form − → A . − → p involving the vector potential and momentum operators, in contrast to traditional implementations for isolated atomic and molecular systems 5,6,8 where the length-gauge form − → r . − → E coupling the position and electric field operators is employed. Numerical implementations of velocity-gauge RT-TDDFT (VG-RT-TDDFT) based on real-space grids 26,27 and full potential linearized augmented planewaves (FP-LAPW) 28 have been demonstrated in recent years and used to investigate laser-induced valence electron-and spin-dynamics 11,13 and related phenomena in solid-state systems with encouraging results.
Valence electron excitations which involve energy scales of a few eV are amenable to a very efficient treatment using uniform real-space based RT-TDDFT 26,27 as the valence electron density and its fluctuations are smoothlyvarying in space and can be well represented on relatively coarse grids with grid spacings on the order of 0.5 a.u. Analogously, in planewave implementations 17 a small kinetic energy cutoff on the order of a few tens of Rydberg is sufficient.
In practical simulations that are primarily concerned with valence electron dynamics, one often employs either a pseudopotential or frozen-core approximation whereby the localized inner-shell core electrons are effectively eliminated from the description facilitating low-cost real-space or planewave expansions. However, if one is interested in higher energy excitations on the scale of a few tens to hundreds of eV, especially in the context of inner-shell spectroscopies 29 , the core-electrons cannot be disregarded and the much denser real-space grids or higher kinetic energy of planewaves necessary to describe the highly-localized inner-shell orbitals in turn significantly increase the computational cost of such simulations. Adaptive non-uniform real-space grids 30 and FP-LAPW 28 methods afford possible means to circumvent this issue but an alternate approach is to employ a localized basis-set framework. A linear-combination of atomicorbitals (LCAO) 31 approach allows naturally for an efficient treatment of electronic states localized near the atomic nucleus but while potentially sacrificing some variational freedom in the regions of low electron density away from nuclei especially when highly diffuse electronic states are involved. Nevertheless, such an approach might represent a worthwhile compromise in the context of theoretical simulations aimed at understanding inner-shell spectroscopies employing extreme ultraviolet (XUV) 32 or X-ray radiation 29 .
The advent of X-ray free electron lasers (FELs) 33 and of high harmonic generation based attosecond XUV laser pulses 32 has led to the development of novel ultrafast spectroscopies that utilize inner-shell excitations to investigate electron dynamics on femtosecond and sub-femtosecond time scales 34,35 . In particular, attosecond XUV spectroscopy has been extended to solid-state systems exploring for instance, the early-time dynamics of electrons excited across the band gap of semiconductor materials irradiated by intense few femtosecond near-infrared (NIR) pulses 36,37 . More recently attosecond time-resolved core-exciton dynamics have also been investigated in a solid-state context 38 .
As experimental capabilities utilizing short XUV or X-ray pulses advance further into areas such as non-linear X-ray spectroscopy, theoretical tools that can treat valence and core electron dynamics efficiently and on the same footing are necessary 39 . All-electron Gaussian-type orbital (GTO) basis-set implementations of RT-TDDFT based on the standard length-gauge have already been demonstrated 8 and have been utilized for simulating inner-shell spectroscopies in molecular systems 40 . However, to our knowledge, such GTO RT-TDDFT implementations have not yet been extended to include periodic boundary conditions and the velocity-gauge formalism for treating condensed-phase systems.
It is in this scenario that we explore a velocity-gauge implementation of RT-TDDFT within an LCAO basis-set framework with the aim of assessing the level of agreement with real-space grid or FP-LAPW methods that such an approach can facilitate. The VG-RT-TDDFT implementation described here is incorpo-

Formalism
The velocity-gauge RT-TDDFT formalism implemented in this work is based on the one due to Bertsch et al 4,11 . In RT-TDDFT, the time-dependent Kohn- whereĤ KS and ψ i ( − → r , t) are the Kohn-Sham (KS) single-particle Hamiltonian and orbitals respectively, are integrated in the time-domain. Using the KS Hamiltonian in the length-gauge normally adopted with finite systems, the TDKS equations for electron dynamics take the form whereV ion is the electron-ion interaction, the integral represents the Hartree yielding the velocity-gauge TDKS equations wherein the vector potential − → A (t) appears in the kinetic term. In the case where non-local pseudopotentials are used, the gauge field also transforms the electron-ion interaction term toV ion related toV ion bŷ Although the length-gauge and velocity-gauge forms yield the same result for finite systems, the Hamiltonian in equation 5 is periodic for spatially uniform external electric fields, which allows a Bloch representation to be used in solving the velocity-gauge TDKS equations. Integrating equation 5 in time yields the time-dependent electron density n( − → r , t) = i |ψ i ( − → r , t)| 2 as well as the time- where the time-dependent current density features the generalized momentum which accounts for the possible use of nonlocal pseudopotentials. Once the timedependent density and current density are available frequency domain quantities can be accessed through Fourier transforms 2,4 . For instance, an impulsive field − → where θ(t) is the Heaviside step function, can be applied and the time-dependent current Fourier transformed to yield in the linear response regime the frequency dependent conductivity and the frequency dependent dielectric function In equation 10, f (t) represents a filtering function inserted to avoid an abrupt cutoff of the integrand at end of the time-evolution period T .

Implementation in SIESTA
The KS Hamiltonian on the right hand side of equation 5 differs from its length-gauge counterpart in the structure of the kinetic and non-local pseudopotential terms. The Hartree and semi-local exchange-correlation (XC) terms are structurally unaffected by the velocity-gauge transformation and within the adiabatic approximation, can be calculated using standard procedures laid out for ground state simulations 41 . Expanding the kinetic term that now includes the time-dependent vector potential yields In constructing the Hamiltonian matrix elements H mn = m|H|n over SIESTA basis functions (m, n), we note that the first term on the right hand side involving ∇ 2 is simply the standard kinetic energy term. The second term in the expression requires evaluating matrix elements of the − → ∇ operator in the LCAO basis. These are also obtained using the standard reciprocal space scheme for two-center integrals in SIESTA 41 . Matrix elements of the third term involv- is spatially homogenous within the present description.
The vector potential entering the modified nonlocal pseudopotential operator in equation 6 gives rise to two-center integrals of the form where − → R mn is the vector connecting the nuclear coordinates of the two atom- so that the time dependence enters via the expansion coefficients c k im . Adopting atomic units, the TDKS equation 5, is then recast in terms of the coefficient where H k and S k represent the Bloch-downfolded Hamiltonian and overlap matrices in the SIESTA basis. Efficient unitary time evolution of the coefficient vectors c k i is carried out using the Crank-Nicholson 45 scheme whereby Additionally, we utilize parallelism over k-points through both distributed and shared memory parallelization in order to speed up the simulations on multiprocessor architectures.

Linear response and laser induced dynamics in bulk Silicon
With the aim of comparing the LCAO basis implementation of VG-RT-TDDFT within SIESTA with existing real-space implementations, we carry out a number of simulations on prototypical bulk Silicon investigating both linear response and strong laser field induced dynamics. As the real-space grid VG-RT-TDDFT implementation, we employ the well-established ARTED code 27 which is a precursor to the SALMON project 46 , and has been used extensively for studying laser induced valence electron dynamics in solid state systems 11,18,22 .  Next we investigate the interaction of intense laser fields with bulk Si. For this purpose, we consider a 10 fs infrared (IR) laser pulse with a carrier frequency of 1.6 eV and a sin 2 envelope centered at t = 5 fs as shown in Fig.2(a).
The electric field of the laser pulse is oriented along the z-axis of the conventional unitcell and we assume transverse boundary conditions so that surface polarization effects do not play a role 11 . The action of this pulse for field intensities ranging from 10 10 -10 13 W/cm 2 (Fig.2(a)) is simulated within both SIESTA and ARTED.
In Fig. 2 The agreement between SIESTA and ARTED in the frequency domain is also satisfactory over the range of intensities considered. We note in particular that high-harmonic signals from the two codes compare well even though we do not employ very diffuse functions within the LCAO basis set. This is because in a solid state system like bulk Si, regions of extremely low electron density such as in the tails of molecular wavefunctions decaying into vacuum do not occur and the multiple-ζ basis is able to adequately describe density fluctuations.
Having analyzed the laser driven current in both the time and frequency domains, we now consider the energetics of laser-matter interaction. As the laser pulse passes through bulk Si, electrons can be excited across the band gap leading to net energy absorption by the material. We note that, in the present instance, the central frequency of the exciting field at ∼ 1.6 eV, is well below the

Core-level spectroscopy in Silicon Carbide
The primary convenience offered by an LCAO based RT-TDDFT approach over real-space grid basis methods is with respect to the treatment of localized core orbitals relevant to extreme ultraviolet (XUV) and X-ray spectroscopies.
In this subsection we therefore investigate core-level response in 2H Silicon Carbide (SiC) (Fig. 4(a)) to demonstrate the utility of the present LCAO approach in this context. In particular we consider the L-edge of Si and the K-edge of C in the XUV and soft X-ray ranges respectively. Semi-core and core states can be treated within SIESTA by constructing appropriately modified pseudopotentials 31 . For instance, in the case of the C K-edge (Si L-edge), we prepare a pseudopotential that explicitly pseudizes the C:{1s, 2p, 3d} (Si:{2s, 2p, 3d}) states. The C 2s (Si 3p) state is then automatically obtained as a higher energy solution of the atomic Schrodinger equation and its wavefunction has the correct nodal structure due to orthogonality with the C 1s (Si 2p) wavefunction (See Fig. 4(b)). The C 1s (Si 2p) wavefunction is naturally consistent with the nodelessness criterion of standard pseudization techniques. Therefore, with an appropriate choice of the pseudopotential matching radius, both the C 1s, 2s (Si 2p, 3p) pseudo wavefunctions and eigenvalues closely resemble the corresponding all-electron counterparts as shown in Fig. 4(b). Similar considerations apply for the 2s and 3s states of Si. Such a pseudopotential can be incorporated into the SIESTA framework to model specific semi-core states in solid-state and molecular simulations.
Within Elk, a 24x24x40 planewave grid in conjunction with 28 local orbitals is utilized as the basis set. In both codes, a 20x20x12 k-point grid is used for Brillouin zone sampling. The electronic density of states (DOS) of SiC calculated from SIESTA and Elk are shown in Fig. 5(a) over an extended energy range. As one would expect, very good agreement is obtained for the DOS in a range spanning several eV on either side of the Fermi level. For the Si 2p (C 1s) core states near -91 eV (-262 eV), we find that KS energy band centers from the two codes agree to within 0.15 eV (0.06 eV) with the SIESTA eigenenergies slightly blue-shifted. We note that this difference on the order of 0.1% is small compared to other sources of uncertainly within the theoretical approach such as due to the XC approximation. Furthermore, since spin-orbit coupling is neglected in the present simulations, no splitting of the Si 2p core-level DOS is seen. In Fig. 5(b) the imaginary part of the frequency dependent dielectric function ( zz (ω)) as obtained from SIESTA and Elk, calculated using the independent-particle approximation (IPA) for light polarization along the c-axis of the crystal (See Fig. 4(a)) is compared over a 280 eV energy range. Spectral features originating from relevant low-energy valence and high-energy core excitations show satisfactory agreement both with respect to frequencies and relative oscillator strengths especially at energies near the absorption edges. With regards to core-level spectra, IPA derived frequencies show small differences on the order of 0.1 eV consistent with those in the underlying KS eigenvalues.
Furthermore, while spectral features are in good agreement overall, absolute oscillator strengths at higher energies are somewhat larger in the LCAO spectra at the Si L-edge. We tentatively assign this to differences in the description of continuum wavefunctions in LCAO and FP-LAPW methods while noting that in the latter approach, oscillator strengths at high energies depend somewhat on the choice of muffin-tin radii. These differences nevertheless should not affect spectral interpretation in the near-edge region.
Next we investigate real-time response in SiC using the VG-RT-TDDFT implementations in SIESTA and Elk with a particular emphasis on core-level excitations at the Si L-and C K-edges. Relevant numerical parameters for realtime simulations are the same as the ones adopted to calculate the DOS and IPA response. Within Elk, a total of 115 KS states, with 18 occupied and 97 empty bands are included in the real-time description. Since core excitations at high energies are characterized by rapid oscillations of the electron density, a smaller time step is typically required to propagate the TDKS equations compared to valence-only simulations. Accordingly, we employ a time step of 0.01 a.u in both codes (See also discussion around Fig. 7). By the same token however, time propagation does not need to be carried out for long periods to sample the fast core-electron oscillations, so the system is propagated in this instance for a total of 10 fs. In Fig. 6, the time-dependent current density J z (t) induced in response to a weak 0.001 a.u impulsive electric field applied at time zero along  Fig. 6(a) shows a magnified view of J z (t) from SIESTA during the interval between 2 -2.5 fs from which it is apparent that the current density is characterized by high frequency core-electron oscillations of a smaller amplitude superimposed imposed upon slower valence oscillations. As before, J z (t) can be Fourier transformed to calculate the zz-component of the linear dielectric function zz (ω) in the frequency domain whose real and imaginary parts are plotted in Fig 6(b,c).
The real part (Fig 6) is plotted over a 20 eV range in the valence region and exhibits good agreement between the LCAO and FP-LAPW approaches. The imaginary part of zz (ω), which is relevant to absorption spectroscopies, is plot- since the valence-excitation background at X-ray energies is not very sensitive to chemical differences and in the context of time-resolved spectroscopies, does not contribute significantly to the differential pump-probe signal.
Finally, we discuss briefly the issue of convergence of the core-level spectra  with integration time step. In Fig. 7 we show the C K-edge spectra obtained from SIESTA based VG-RT-TDDFT for three different integration time steps between 0.01 -0.04 a.u. At all three time steps, the Crank-Nicholson integration scheme is stable over the 10 fs propagation time. At the smallest time step of 0.01 a.u the absolute energy positions of the spectral features are converged to within 0.1 eV suggesting that the same is adequate to sample rapid oscillations due to C K-edge derived excitations. In contrast the spectra for time steps of 0.02 a.u and 0.04 a.u are shifted to higher energies. Nevertheless, it is apparent that while the lower sampling rate affects absolute energy positions, the relative energy differences between different spectral features are practically identical for all three sampling rates. This provides a numerical illustration of the approximate separability of the electron dynamics into fast core and slow valence oscillations with the former primarily determining the absolute energy positions of the core-excitation edges and the latter being more relevant for the detailed structure of the near-edge features. Therefore, in some applications where reproducing the absolute core-excitation edge energy is not crucial, larger time steps can in principle be used to simulate near-edge spectral features. A similar separation of time scales argument has also been utilized by Lee et al for the simulation of core-level spectra within the local time-correlation approach 50 .

Conclusions and Outlook
In summary, we have implemented the velocity-gauge formalism of realtime TDDFT within a numerical atomic orbital based first-principles framework and carried out a comparative assessment of time-varying electric field response properties against results from real-space grid (RSG) and FP-LAPW codes. The overall agreement between the current LCAO basis approach and RSG or FP-LAPW methods is satisfactory in that the basis set or pseudopotential driven differences in the observable quantities simulated are small compared to errors intrinsic to the underlying first-principles density functional approximations. LCAO and RSG/FP-LAPW predicted excitation frequencies differ negligibly in the valence region and by ∼0.1 eV in the core-region. The energetics of laser-matter interaction are also well described by the LCAO framework and are within 5% of RSG results. The fact that this level of agreement is achieved without employing very large basis sets or targeted optimization of the basis functions suggests that there is potentially scope for obtaining even closer agreement where necessary, with some additional effort. The primary advantages afforded by the LCAO approach are the relatively small size and shortranged character of the basis set which leads to efficient representations for the density and Hamiltonian matrices and favorbale scaling with increasing system size 41 . The present LCAO implementation of VG-RT-TDDFT in-principle retains these advanages while enabling the study of time-resolved and nonlinear XUV/X-ray spectroscopies at similar computational cost as simulating valenceonly dynamics. However, in the context of modeling laser-matter interaction in small unitcells with dense reciprocal space-sampling while employing adiabatic semi-local XC functionals, it does not afford significant advantages over more mature existing FP-LAPW based VG-RT-TDDFT 13 . While some efficiencies could potentially be harnessed in supercell simulations of low dimensional systems that include large vacuum regions, in order to fully take advantage of the small basis set sizes inherent to the LCAO approach, the implementation must be extended to make use of non-local XC functionals 51 to describe excitonic effects and density-matrix evolution to incorporate coupling to external baths within a Liouville picture 2,52 . Efforts along these lines are currently underway.