Semi-classical description of electron dynamics in extended systems under intense laser fields

We propose a semi-classical approach based on the Vlasov equation to describe the time-dependent electronic dynamics in a bulk simple metal under an ultrashort intense laser pulse. We include in the effective potential not only the ionic Coulomb potential and mean-field electronic Coulomb potential from the one-body electron distribution but also the exchange-correlation potential within the local density approximation (LDA). The initial ground state is obtained by the Thomas-Fermi model. To numerically solve the Vlasov equation, we extend the pseudo-particle method, previously used for nuclei and atomic clusters, to solids, taking the periodic boundary condition into account. We apply the present implementation to a bulk aluminum (FCC) conventional unit cell irradiated with a short laser pulse. The optical conductivity, refractive index, extinction coefficient, and reflectivity as well as energy absorption calculated with the Vlasov-LDA method are in excellent agreement with the results by the time-dependent density functional theory and experimental references.

We propose a semi-classical approach based on the Vlasov equation to describe the time-dependent electronic dynamics in a bulk simple metal under an ultrashort intense laser pulse. We include in the effective potential not only the ionic Coulomb potential and mean-field electronic Coulomb potential from the one-body electron distribution but also the exchange-correlation potential within the local density approximation (LDA). The initial ground state is obtained by the Thomas-Fermi model. To numerically solve the Vlasov equation, we extend the pseudo-particle method, previously used for nuclei and atomic clusters, to solids, taking the periodic boundary condition into account. We apply the present implementation to a bulk aluminum (FCC) conventional unit cell irradiated with a short laser pulse. The optical conductivity, refractive index, extinction coefficient, and reflectivity as well as energy absorption calculated with the Vlasov-LDA method are in excellent agreement with the results by the time-dependent density functional theory and experimental references. The interaction of ultrashort (fs-ps) intense laser pulses with solids is relevant to a wide area of research ranging from high-harmonic generation [1][2][3][4][5] to material machining [6][7][8][9][10][11][12][13][14]. The process of ultrafast laser micromachining, which can suppress heat-affected zone, starts from the energy transfer from the laser to the material by electron excitation, followed by that from the hot electrons to the lattice. As a result, the material undergoes phase and/or structural transition [15], leaving a change of the optical constants or a defect behind [16], which eventually leads to ablation, drilling, or structuring [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].
The comprehensive modeling of laser material machining is highly complex, multi-scale in both time and space, multi-phase (solid, fluid, plasma, cluster, etc.), and possibly accompanied by chemical reactions. Plasma or continuum models [6,8,[33][34][35][36], for example, have been employed to describe and simulate such processes, advancing understanding. However, they have difficulties in examining initial transient dynamics before the local thermodynamic equilibrium is reached.
It has now become possible to describe the attosecondfemtosecond electron dynamics under intense laser fields with the time-dependent density-functional the-In TDDFT, each electron orbital satisfies the timedependent Kohn-Sham (TDKS) equation [see Eq. (1) below]. The leading order of a semiclassical expansion of the TDKS equation reduces to the Vlasov equation, which describes the temporal evolution of the electron distribution function in phase space. Thus, Vlasov-based approaches are expected to be a cost-effective alternative to TDDFT, in particular, for metals. Such approaches have previously been applied to ionization and explosion dynamics of molecules [44,45] and metal clusters [46][47][48][49][50][51]. The Vlasov equation is numerically solved with so-called pseudo-particle methods in these studies, which represent the electron cloud as an assembly of classical test particles whose motion is governed by Newton's equations of motion. There are several reports of application to Na clusters, well agreeing with TDDFT results [49][50][51][52][53].
In this paper, we extend the pseudo-particle method based on the Vlasov equation to the description of electron dynamics in extended systems under intense laser fields. The effective potential acting on the electrons contains not only the ionic potential, interelectronic Hartree potential, and interaction with laser but also the exchange-correlation potential within the local-density approximation (LDA), and incorporates the periodic boundary condition. We apply the present method to bulk aluminum. The calculated optical conductivity, refractive index, extinction coefficient, and reflectivity as well as energy absorption are in excellent agreement with TDDFT calculations and experimental references.
The present paper is organized as follows. Section II describes our simulation methods. We review the Vlasov equation and describe our numerical implementations with the periodic boundary condition. In Sec. III we describe numerical application to bulk aluminum and compare the results with TDDFT and measurement values. The conclusions are given in Sec. IV.

A. Vlasov equation
Among the methods for treating quantum many-body dynamics, TDDFT provides a feasible computational framework for treating electronic systems' optical response or charged particles' collision phenomena [54]. The time propagation of a N e -electron system comes down to solving a set of equations for the Kohn-Sham orbitals {φ i (r, t)} that evolve in a self-consistent mean field [55], where, denotes the Kohn-Sham Hamiltonian, m the electron mass, V eff the effective potential (see below), and the time-dependent electron density n e (r, t) is defined as, Analogously, the density operatorρ(t) is defined as, whose evolution is govenrned by the von-Neumann equation (vNE), Performing the Wigner transformation [56] and taking the limit → 0, the density operatorρ(t) is mapped onto a real function f (r, p, t), which obeys the Vlasov equation, which is a classical alternative to the vNE, Eq. (5), where p is the electron canonical momentum. Here, f (r, p, t) is interpreted as the electron distribution in phase space. The effective potential V eff is a functional of the electron density distribution n e (r, t) and decomposed into, with the exchange-correlation potential V xc , external field potential V ext , and, (8) where i, V ps and V H denote the label of ions and the spherically symmetric ionic pseudopotential and the electron-electron Hartree potential, respectively.
Several previous works for Na clusters have used their original pseudo potentials [46,47], adjusted so that the simulation results reproduce the static and dynamical properties of the system. In this work, instead, we employ the modified Heine-Abarenkov type local pseudo potential for V ps , where z is the number of the valence electrons, and A, R, α, and β are material dependent parameters determined by ab initio density functional formalism in Ref. [57], thus, independently from Vlasov simulations. Their values for the bulk aluminum crystal are A = 3.574 a.u., α = 3.635 a.u., β = 0.8343 a.u., R = 0.334 a.u., z = 3. V H is evaluated by solving the Poisson equation, Here, let us introduce a real-space simulation box Ω, on which the periodic boundary condition is imposed, and translation vectors G. Ω is defined as, where {e j } are the lattice vectors along the j-axis (j = x, y, z), whose lengths are denoted by L j = |e j |. Integrals with respect to r are taken over Ω in what follows. The translation vectors are given by, Taking the periodic boundary condition into account, the Coulomb terms V ps and V H are represented as a Fourier series expansion. The pseudo potential term is rewritten as, where N ion denotes the number of ions in Ω and, with Q being the coordinates in the Fourier domain (Q = |Q|), F [·] the Fourier series expansion within Ω, and, One obtains the solution of the Poisson equation (Eq. (10)) for the electron density n e given within Ω as, For the exchange-correlation potential V xc , one employs the LDA by Perdew and Zunger [58]. The laserelectron interaction is described in the length gauge, within the dipole approximation, where E denotes the laser electric field vector. In this case, p becomes the kinetic momentum, and, thus, the electronic current density J(t) averaged over Ω is given by, B. Numerical implementations

Pseudo-particle method
The direct propagation of the distribution function would require the treatment of six-dimensional timedependent function on grids [59]. To avoid such a massive computation, one introduces the pseudo-particle method [46][47][48]60], where the distribution function f (r, p, t) is expressed by a set of pseudo-particles with mass m as, Here r i , p i are the position and canonical momentum of each pseudo particle labeled by i. The total number of pseudo particles N pp is given by N pp = N s N e , where N s and N e are the number of pseudo particles per electron and the total number of the electrons contained in Ω, respectively. Statistical error is reduced by increasing N s . N s is set to 10000 in this study. g r (r) and g p (p) denote smoothing kernel functions for the position and momentum, respectively, of Gaussian forms, where d r and d p are smoothing widths. Only the nearest neighbor cells are included in summation over G in Eq. (20). The kernel functions are normalized as, g p (p)dp = 1, so that, The scattering cross-section of the electron and the effective potential is adjusted through d r ; the smaller d r , the larger the cross-section. Here we use d r = 0.575 a.u. so that the time-dependent energy absorption reproduces the TDLDA results. In the present collision-less case, d p is not used explicitly.
The field quantities such as V eff and n e are evaluated on three-dimensional grids discretized into N j (j = x, y, z) intervals on the j axis with a spatial step ∆j = L j /N j . In our calculation we set ∆x = ∆y = ∆z = 0.5 a.u. Here, d r /∆x ≃ 1.15 is a good parameterization leading to stable simulation [47]. It should be noticed that d r is the only adjustable parameter in our formalism. The electron density on a grid point r is calculated as, The current density J(t) [Eq. (18)] is evaluated as, The Hamiltonian in pseudo-particle representation is written as, The motion of each pseudo particle is governed by the Newton equations under the effective potential V eff with the periodic boundary condition as, As long as pseudo-particle canonical variables r i , p i obey the Newton equation (Eq. (28)), one-body distribution (Eq. (19)) satisfies the Vlasov equation (Eq. (6)). The force term is given as the gradient of the N pp -body Hamiltonian H pp . One numerically integrates it as, using the analytical form of ∇ ri g r (r i − r), The integration of Eq. (28) is performed by the Verlet method [61] with time step ∆t = 0.02 a.u. Particles exiting Ω are to re-enter Ω from the other side.

Ground state
The initial state is the stationary solution of the Vlasov equation is variationally minimized with respect to n e (r) under the constraint that the box Ω contains N e electrons. This leads to the following coupled equations, where µ denotes the chemical potential, playing the role of a Lagrange multiplier. These equations are to be solved for n e (r) self-consistently. An adopted algorithm to solve the coupled equations Eq. (32) and Eq. (33) is shown in Fig. 1.
First, the chemical potential µ and the electron density n in e in the real space are guessed so that the total 1: procedure Ground state preparation 2: (-Initialization-) 3: initial guess of µ and n in e 4: 5: (-Self-consistent determination of µ and ne-) 6: while ∆n > ǫ (ǫ = 10 −7 ) do 7: set pseudo-particle position {ri} 8: {ri} → nps(r) using Eq. (25) 9: nps → V eff [nps](r) using Eq. (33) 10: V eff [nps](r) → ne 11: find appropriate µ using random numbers (line 7, also see below). The electron density n ps realized by the pseudo-particle distribution is calculated through Eq. (25) (line 8). The effective potential V eff (r) is obtained by substituting n ps into the right-hand side of Eq. (33) (line 9). Then, we update the electron density n e (r) by substituting thus obtained V eff (r) to the Eq. (32) and solving it with respect to n e (line 10), simultaneously updating µ by a bisection method to satisfy the condition that the total number of electrons is N e (line 11). The updated n e is used as n in e in the next iteration of the loop (line 13). One repeats the above operations till convergence, Ω dr|n in e −n e | < ǫ, where we set ǫ = 10 −7 here for crystalline Al (line 12). After convergence, one distributes the momenta of the pseudo particles uniformly within the local Fermi radius p f by acceptance-rejection sampling of uniform pseudorandom numbers (lines [17][18][19][20][21][22][23][24][25]. The algorithm to distribute the pseudo particles (line 7 in Fig. 1) is shown in Fig. 2.
We introduce sub-grid points (line 3) by dividing each voxel of the computational grid into N j inp regions along the j-axis (j = x, y, z). The electron density n inp e on each sub-grid point is evaluated, based on the trilinear interpolation from those of the surrounding eight computational grid points, from which one calculates the number of pseudo particles N local pp around the sub-grid point (lines 4-5). Then, the N local pp pseudo particles are uniformly distributed around the sub-grid point using random num-1: procedure how to set pseudo particle position

C. Linear response
We evaluate the linear optical response via impulse response by doing dynamical simulations with the initial pseudo-particle momenta p i shifted from the groundstate values p GS i by a small amount ∆p, where ∆p = (0, 0, 0.1 a.u.) in this study. This is equivalent to the application of an impulse electric field, where δ(t) is the delta function. Noting that this field has a constant power spectrum across all frequencies, one can readily obtain the optical conductivity as, where ∆p m andĴ m (m, n = x, y, z) denote the m component of the momentum shift and the temporal Fourier transform of the current density, respectively. The fast Fourier transformation algorithm [62] is used for the evaluation ofĴ m (ω). Assuming isotropic media, the dielectric function ε mm (ω), the complex refractive index n(ω), and the reflectivity R(ω) are given by, respectively, especially, ε xx (ω) = ε yy (ω) = ε zz (ω).

III. RESULTS
In this section, we compare the results of the Vlasov-LDA simulations for extended systems described in the previous section with the experimental values and the TDDFT results obtained by the open source code SALMON [37,63,64]. We take aluminum as a target material. For Vlasov-LDA, simulation parameters are N s = 10000, N e = 12, and time step ∆t = 0.025 a.u.. For TDDFT, we employ a norm-conserving pseudopotential [65] and the LDA functional [58], with the number of k-points 48 3 , number of real-space grids 14 3 , and dt = 0.15 a.u. We assume an external electric field linearly polarized along the Γ − X direction of the following temporal profile: where E 0 denotes the field amplitude, ω the photon energy, and T the (foot-to-foot) full pulse duration. The corresponding full width at half maximum duration of the laser intensity profile is about 0.36T .

A. Linear response
Let us first discuss the complex optical conductivity, refractive index, extinction coefficient, and reflectivity as a function of photon energy. Despite the simpleness of the Vlasov-LDA approach, its results excellently agree with the TDDFT results and experimental values (Fig. 3), especially above 2 eV photon energy. The peak and dip around 1.5 eV in the TDDFT results are due to interband absorption, which is not reproduced by the present Vlasov approach, since the latter takes only the single free-electron dispersion into account. Focusing on reflectivity behavior around the plasma frequency, one finds some differences between the two approaches. This difference would be contributions by the abovementioned interband resonance and non-unity effective mass in TDDFT. We have confirmed it through the decomposition of the response obtained by TDDFT into Drude and Lorentz model components. With the resonance energy set to 1.85 eV, the biggest oscillator around 1.5 eV [67], The estimated effective mass m eff is 1.09m, and the damping constant is 0.51 eV −1 , consistent with the values reported previously (1.16m [68] and 0.80 eV −1 [67], respectively). The loss functions, Imǫ(ω) −1 , are shown in Fig. 4. The TDDFT result is excellently reproduced by the combined Drude and Lorentz contributions with m eff = 1.09m. Although Vlasov-LDA overestimates the plasma frequency compared to TDDFT, it agrees with the Drude model with m eff = m.

B. Energy absorption
Let us next investigate the energy absorption from the laser pulse. We evaluate the energy absorption by the electrons as their energy increment by the pulse irradiation. The energy is calculated as ∆E = H pp (t = ∞) − H pp (t = 0) in the Vlasov-LDA simulation and as h KS (t = ∞) − h KS (t = 0) in the TDDFT simulation. We show the fluence dependence for the fixed intensity (10 12 W/cm 2 ) at 80 nm wavelength in Fig. 5 and that for the fixed pulse width of 3.8 fs at 200 and 400 nm wavelengths in Fig. 7.
We see in Fig  in the lower fluence region ( 50 mJ/cm 2 ). On the other hand, the Vlasov-LDA does not reproduce the TDDFT results for the higher fluence, where the latter deviate from the linear behavior and even decrease with increasing fluence. This difference is due to Rabi-like oscillation [69], as confirmed in Fig. 6, which shows the temporal evolution of absorbed energy for several pulse widths. The maximum electron energy gain during the pulse, which are indicated by horizontal dashed lines, does not depend much on the pulse width, suggesting Rabi-like coherent oscillation. Thus, there is an optimum pulse width for a fixed intensity in terms of energy absorption. Figure 7 indicates that the energy absorption calculated by the Vlasov-LDA approach exhibit a linear dependence on fluence or pulse intensity in the whole range. We can see a nonlinear behavior, on the other hand, in the TDDFT results in the higher fluence range ( 40 mJ/cm 2 for 200 nm and 2 mJ/cm 2 for 400 nm. This would be interpreted as saturable absorption as is widely observed in various materials [70]. We could not obtain the Vlasov-LDA results for the low fluence region ( 2 mJ/cm 2 ) because of statistical error. This could be improved by increasing the total number of pseudo particles, in principle. Nevertheless, the Vlasov-LDA results, if extrapolated to the low fluence region, appear to agree well with the TDDFT results. Figure 8 shows the temporal evolution of the current density and absorbed energy for 10 12 W/cm 2 peak intensity, 200 nm wavelength, and 3.8 fs pulse width. Again, overall, the Vlasov results excellently reproduce the TDDFT results. In Fig. 8 (b), although energy fluctuation due to the pseudo-particle statistical error is seen  at < 2 fs, it becomes negligible at the end of the pulse.
Our code is partially parallelized using OpenMP and MPI. One of the most time-consuming parts is the Fourier transformation, which is computed by naive approach. Nevertheless, the computational time of the present Vlasov-LDA code is typically only 1/20 of that of TDDFT using the SALMON code. With more sophistication and parallelization, the efficiency of the Vlasov-LDA method will be further improved, which will be advantageous for applications such as parameter optimization in laser material processing.

IV. CONCLUSIONS
We have extended the Vlasov-LDA semi-classical approach and implemented it with the pseudo-particle method to periodic systems in order to compute the electron dynamics in solids, especially in metals, under ultrashort intense laser pulses. The Vlasov equation can be regarded as the leading order of a semiclassical expansion of the time-dependent Kohn-Sham equations. The electronic distribution function is expressed by pseudoparticles, incorporating the periodic boundary condition. They play the role of Lagrangian markers embedded randomly in the electron gas. The initial distribution is calculated from the Thomas-Fermi model.
We have applied this approach to crystalline aluminum. Although the method has only one adjustable parameter d r , the calculated optical conductivity, refractive index, extinction coefficient, and reflectivity as well as energy absorption are overall in excellent agreement with the TDDFT and experimental results over a wide range of photon energy and fluence, demonstrating the capability of the present approach to accurately describe the dynamics of metallic conduction-band electrons. On the other hand, the Vlasov results deviate from the TDDFT ones around 1.5 eV photon energy, where interband transitions are involved, and at the high-fluence region, where a Rabi-like oscillation takes place.
The next step will be to incorporate electron-electron collisions, the description of which is limited in TDDFT. Vlasov-LDA is expected to provide valuable insights into complex laser-material processing if we further couple it with molecular dynamics, electromagnetic field analysis, and other continuum models.