Ab initio simulation of attosecond transient absorption spectroscopy in two-dimensional materials

We extend the first-principles analysis of attosecond transient absorption spectroscopy to two-dimensional materials. As an example of two-dimensional materials, we apply the analysis to monolayer hexagonal boron nitride (h-BN) and compute its transient optical properties under intense few-cycle infrared laser pulses. Nonadiabatic features are observed in the computed transient absorption spectra. To elucidate the microscopic origin of these features, we analyze the electronic structure of h-BN with density functional theory and investigate the dynamics of specific energy bands with a simple two-band model. Finally, we find that laser-induced intraband transitions play a significant role in the transient absorption even for the two-dimensional material and that the nonadiabatic features are induced by the dynamical Franz-Keldysh effect with an anomalous band dispersion.


Introduction
Attosecond transient absorption spectroscopy (ATAS) is a technique that employs attosecond laser pulses to probe modifications to the optical absorption in the time domain. Because its attosecond time-resolution is shorter than a typical time-scale of electron dynamics, it can naturally capture ultrafast electron dynamics in matter. Therefore, ATAS is one of the key experimental techniques to explore ultrafast phenomena where nonlinear electron dynamics plays a significant role. Indeed, ATAS has been extensively applied to atoms and molecules over the past decade to investigate ultrafast electron dynamics in such relatively small systems [1][2][3][4][5][6]. Recently, ATAS has been further extended to solid-state materials to investigate rather complex electron dynamics in various phenomena such as band-gap renormalization [7], petahertz optical drive [8], dynamical Franz-Keldysh effect [9], photocarrier relaxation [10], exciton dynamics [11] as well as photocarrier injection [12].
Two-dimensional (2D) materials have been attracting great interest from both fundamental and technological points of view because of their remarkable properties. For example, graphene-like 2D materials such as hexagonal boron-nitride (h-BN) have been intensively studied for applications in optelectronic devices and energy storage among others [13][14][15]. Transition metal dichalcogenides (TMDs), another kind of 2D materials, are also attracting attention because of their interesting properties deriving from strong spin-orbit coupling and a large direct-bandgap [15][16][17]. Thanks to recent development of laser technologies, laser-induced ultrafast electron dynamics in such 2D materials has been intensively studied both experimentally and theoretically [18][19][20][21][22][23][24]. For example, the high-order harmonic generation was recently studied in 2D materials [25][26][27][28]. Surprisingly, it was demonstrated that elliptically polarized light can enhance the high-order harmonic generation in graphene [25].
In order to fully understand such laser-induced ultrafast electron dynamics, ATAS is one of the most desirable experimental techniques as it can resolve the dynamics with natural time-resolution. While ATAS for solid-state materials provides a wealth of information on microscopic electron dynamics in materials, the resulting experimental data are often very complicated due to complex electronic structure and highly-nonlinear effects in the dynamics. Therefore, it is hard to directly interpret the experimental results. Ab-initio simulation based on the time-dependent density functional theory (TDDFT) [29] is a powerful tool to describe such complex electron dynamics and to provide microscopic insight into the phenomena. Indeed, ab-initio TDDFT simulations have been applied to ATAS experiments for solid-state systems and played a significant role to construct the interpretation [7,9,10,12].
Likewise, experimental results of ATAS for 2D materials are expected to be rather complicated. Therefore, ab-initio TDDFT simulations are a good candidate to analyze ATAS for 2D materials. To realize such simulation methods, we extend the TDDFT pump-probe simulation [30] of transient absorption spectroscopy to 2D materials. To demonstrate the ab-initio ATAS simulation, we employ monolayer h-BN as an example of a 2D material and investigate its transient optical property under intense infrared (IR) laser pulses. Furthermore, we analyze the obtained transient absorption spectra of h-BN with a 2D parabolic two-band model and clarify the microscopic mechanism of the ultrafast modification of the optical property.
This paper is organized as follows: In Sec. 2 we first describe the electron dynamics simulation based on TDDFT. Then, we further describe the linear response calculation and the pump-probe simulation to compute static and transient absorption properties, respectively. In Sec. 3 we apply the TDDFT pump-probe simulation to monolayer h-BN and investigate its transient optical properties. Furthermore, we analyze the obtained transient absorption spectra with a parabolic two-band model. Finally, our findings are summarized in Sec. 4.

Methods
In this section, we describe our theoretical and numerical methods to simulate ATAS for 2D materials. As an example of 2D materials, we choose monolayer h-BN in this work.

Electron dynamics simulation for periodic systems
First, we briefly describe the first-principles electron dynamics simulation based on TDDFT [29]. Details of the simulation are described elsewhere [31]. In the framework of TDDFT, electron dynamics in a periodic system is described by the following time-dependent Kohn-Sham equation, where u b k ( r, t) is the Kohn-Sham Bloch wavefunction, b is a band index, and k is the Bloch wave vector. The electron-ion interaction is described by the ionic potentialv ion , while the electron-electron interaction is described by the combination of the Hartree potential v H ( r, t) and the exchange-correlation potential v XC ( r, t). The Hartree potential satisfies the Poisson equation, where the electron density ρ( r, t) is evaluated as with occupation factors, n b k . The exchange-correlation potential v XC ( r, t) is a functional of the electron density. Because its exact form is unknown, one has to approximate it. In this work, we employ the adiabatic local-density approximation [32] (ALDA), where the exchange-correlation potential depends on the electron density locally in space and time. External laser fields are described by spatially-uniform time-varying vector potentials A(t) in the dipole approximation.
To describe monolayer h-BN, we employ three-dimensional periodic boundary conditions with a hexagonal lattice. The real-space atomic configuration in the plane is shown in Fig. 1 (a). In this work, y-axis is taken to be parallel to a B-N bond, while x-axis is taken to be perpendicular to it. The lattice constant of the hexagonal lattice is set to a = 4.76 a.u. The out-of-plane direction is assigned to the z-axis, and the period of that direction is set to L = 50 a.u., which is large enough to avoid spurious inter-layer interactions.
We discretize the hexagonal lattice into three-dimensional real-space grid points with a grid spacing about 0.36 a.u. We further discretize the two-dimensional first Brillouin zone of h-BN into 32 × 32 k-points.
In this work, dynamics of core electrons in the 1s shell of Boron and Nitrogen atoms are not expected to be significant for transient absorption. Therefore, we treat explicitly only valence electrons while the core electrons are frozen. For this purpose, we employ the Hartwigsen-Goedecker-Hutter pseudopotentials [33] to describe interactions between valence electrons and ionic cores.
All the density functional theory (DFT) and TDDFT calculations in this work are carried out with the Octopus code [34].

Optical property from linear response calculation
Here, we revisit the linear response calculation in time domain to obtain optical properties of 2D materials. For this purpose, we consider electron dynamics in monolayer h-BN under an impulsive perturbation, E(t) = e x k 0 δ(t). The polarization of the perturbation is set along the x-direction in Fig. 1 (a), which is an in-plane direction. The strength of the perturbation k 0 is set to 0.005 a.u., which is weak enough to only create excitations on the linear response regime.
In this context, one of the most useful results of TDDFT calculations for 2D materials is the surface density of the electric current, where S is the area of the h-BN sheet in the simulation cell, and v k is the velocity operator: Note that the velocity operator contains a contribution from the ionic potentialv ion because of the nonlocality of the pseudopotential, which is not commutative with the position operator r [31]. Figure 1 (b) shows the surface current density, J(t), for the x-direction as a function of time. Since the impulsive distortion is applied at t = 0, the current is suddenly induced at that time. Then, the induced current shows oscillatory behavior with damping.
One can extract the optical conductivity from the computed current in time domain. The diagonal element of the conductivity can be evaluated as where J(t) and E(t) are the electric current and the electric field for the polarization direction of the impulsive distortion, respectively. To reduce numerical error due to the finite time propagation in the Fourier transform of Eq. (5), we employed a damping factor γ of 0.5 eV/h. Figure 1 (c) shows the real-part of the optical conductivity of h-BN for the in-plane direction, which is computed from the current in Fig. 1 (b). The computed conductivity shows two peaks at around 5 eV and 15 eV. This double-peak structure is consistent with the previous theoretical work based on the first-principles simulation [35]. Because the real-part of the conductivity is closely related to the optical absorption, we treat it as a direct measure of the photo-absorption in this work. We note that the present calculation with ALDA does not properly describe the excitonic effect which is significant in the optical absorption of h-BN around the optical gap [36,37]. However, as will be shown later, we focus on transitions to conduction states with relatively high energy, where the excitonic effect is expected to be less important. Therefore, we simply ignore the excitonic contribution in this work.

Transient optical properties with pump-probe simulations
Now we describe the pump-probe simulation to investigate transient optical properties. Details of the numerical pump-probe simulation are described elsewhere [30].
To compute transient optical properties of materials, we consider electron dynamics under pump E pump (t) and probe E probe (t) electric fields. Solving the time-dependent Kohn-Sham equation (1) under the presence of both the pump and probe fields, one can simulate the electron dynamics induced by the fields. Furthermore, one can calculate the induced electric current by using Eq. (3); we shall call the current under both the pump and probe pulses pump-probe current, J pump−probe (t). Additionally, one can simulate electron dynamics under only the pump pulse, and compute the current: we shall call the current induced solely by the pump pulse pump current, J pump (t). To extract the current induced by the probe pulse under the presence of the pump, we define the probe current, J probe (t), as Finally, referring to the above linear response equation, the transient conductivity can be evaluated as where T probe is the central time of the probe pulse. In contrast to the linear response calculation in equilibrium, the transient conductivity under the pump pulse depends also on the central time of the probe pulse since the time translation symmetry is broken by the pump. Thus, one may investigate dynamics of transient optical properties in time domain by changing the time delay between pump and probe pulses.
To practically carry out the pump-probe simulations, we employ a pump vector potential of the following form; in the domain (−T pump /2 < t < T pump /2) and zero outside of it. Here, E pump is the peak electric field, ω pump is the mean frequency, and T pump is the full duration of the pump pulse.
To perform the pump-probe simulation, we also employ a similar form for the probe vector potential; in the domain (−T probe /2 < t − T delay < T probe /2) and zero outside of it. Here, E probe is the peak electric field, ω probe is the mean frequency, T probe is the full duration of the probe pulse. The time delay between pump and probe pulses is expressed by T delay .

Attosecond transient absorption of monolayer h-BN
In this section, we investigate transient optical properties of monolayer h-BN under intense IR pulses. For this purpose, we employ the TDDFT pump-probe approach explained in Sec. 2.
For the pump-probe simulation, we first set the polarization direction of both the pump and probe pulses to the x-direction in Fig. 1 (a). Thus, the polarization direction of the pump and probe pulses is perpendicular to a B-N bond. Note that h-BN has inversion symmetry along this direction. For the pump pulse, the peak field strength | E pump | is set to 8.7 × 10 8 V/m, the mean frequency ω pump is set to 1.55 eV/h, and the full pulse duration T pump is set to 20 fs. The corresponding full width at half maximum (FWHM) of the pump pulse is about 7.3 fs. For the probe pulse, | E probe | is set to 8.7 × 10 7 V/m, ω probe is set to 15 eV/h, and T probe is set to 1 fs. The corresponding FWHM of the probe pulse is about 260 as. Under these conditions, we repeat the TDDFT pump-probe simulations by changing the pump-probe time-delay, T delay . Figure 2 (a) shows the real-part of the transient conductivity σ T (ω, T delay ) computed by the TDDFT pump-probe simulations. The time-profile of the applied pump electric field is also shown in Fig. 2 (b). In Fig. 2 (a), one sees an oscillatory feature with a V-shaped energy dispersion around 15 eV, which we shall call fishbone structure. As seen from the figure, the frequency of the oscillation in time domain is twice of the pump frequency ω pump . The 2ω pump -oscillation is a direct consequence of the inversion symmetry of the material in the pump-probe direction, which forbids even-order nonlinear responses, because materials with inversion symmetry must show the same response regardless of the sign of the electric field along a symmetry direction. Similar features have been observed in attosecond transient absorption spectroscopy for semiconductors [9,12], and transient terahertz absorption spectroscopy for GaAs quantum wells [38]. These features of bulk materials have been discussed on the basis of the pump-induced intraband transitions and the dynamical Franz-Keldysh effect [39,40]. The dynamical Franz-Keldysh effect is a modification of optical properties of solids under the influence of oscillatory electric fields. It originates from the field-induced intraband transitions via the modification of the dynamical phase factor.  Fig. 1 (a). (b) The time-profile of the applied pump electric field. The result is computed by TDDFT with ALDA.
In order to understand the microscopic origin of the fishbone structure in Fig. 2, we investigate energies of photo-excited electron-hole pairs based on the independent particle approximation. For this purpose, we compute the vertical excitation energies from the highest valence band to the conduction bands. Note that in the independent particle approximation the transition energy from one state to another is identical to the difference of the single-particle energies of those states, c k − v k . Figure 3 shows the computed vertical transition energies from the highest valence band to the conduction bands at each k point along high symmetry paths through Brillouin zone. Transitions around the M point are described in Fig. 3 (a), while those around the Γ point are described in Fig. 3 (b). In Fig. 3 (a), one can see a negative electron-hole mass band around 14 eV at the M point. The negative-mass band is fitted by a parabola (blue-dashed line) to extract band parameters. The extracted reduced electron-hole mass at the M-point is −3.3m e . Furthermore, we numerically confirmed that all the other vertical transitions in Fig. 3 (a) at the M point are dipole forbidden. Therefore, the negative-mass band is expected to have a dominant contribution to the transient absorption spectrum in Fig. 2. To elucidate the contribution from the negative-mass band at the M point, we compute the transient absorption spectrum with a parabolic two-band model. The parabolic two-band model is described in details elsewhere [41].
In the parabolic two-band model, the electron dynamics at each k point in the Brillouin zone is described by the following Schrödinger equation for two-dimensional state vectors, where the off-diagonal matrix element of the Hamiltonian is given by where ∆ cv, k is the transition energy between the valence and the conduction states at k, and p vc, k is the transition momentum. In the Hamiltonian matrix, the Bloch wavevector is shifted based on the acceleration theorem, K(t) = k + e A pump (t)/hc, by the pump pulse [42].
Here, in order to elucidate the role of the pump-induced intraband transitions, only the probe electric field E probe (t) is coupled to the transition momentum p vc , while only the pump vector potential A pump (t) is taken into account in the acceleration theorem. Thus, the probe pulse induces only the interband transitions, while the pump pulse is responsible for intraband transitions.
For this model, we approximate the transition energy by the parabolic band as where g is the band gap, and µ eh is the reduced mass of an electron-hole pair. To describe the negative-mass band at the M point, the reduced mass µ eh is set to µ M = −3.3m e , and the band gap g is set to g,M = 14.1 eV. These parameters are extracted from the parabolic fit shown as the blue-dashed line in in Fig. 3.
The negative-mass band causes an artificial crossing of the conduction and the valence bands. To eliminate the effect from the artificial band crossing of the model, we truncate the transition momentum matrix as follows Here, the transition momentum at the gap p vc is set to 0.1 a.u. Note that the value of the transition momentum does not affect the structure of absorption spectra as long as the field strength of the probe pulse is weak enough to satisfy the linear response condition. We apply a pump-probe simulation to the parabolic two-band model and compute the transient absorption spectrum of the negative-mass band. Since the parabolic-band approximation is an appropriate approximation of the electronic structure only around the band edge, we truncate the signal to reduce the contribution from inappropriate regions by applying the following attenuation function to the transient conductivity Here, g,M is the corresponding band gap at the M-point, which is 14.1 eV. We set the width of the attenuation σ to 1 eV. Figure 4 (a) shows the truncated transient conductivity, σ T (ω, T probe )w M (ω), of the negative-mass band at the M point. One may see that the transient absorption spectrum of the negative-mass parabolic band model in Fig. 4 (a) nicely reproduces the upper-half part of the transient absorption spectrum of the full ab-initio simulation in Fig. 2. This fact indicates that the upper-half part of the fishbone structure originates from this specific band around the M point, which has the anomalous band dispersion with the negative reduced electron-hole mass.  To reveal the other half contribution in the transient absorption spectrum of h-BN, we investigate vertical transitions around the Γ point. Figure 3 (b) shows the vertical transition energies from the highest valence band to conduction bands around the Γ point. Among those of transitions around the energy of 13 eV at the Γ-point, we pick up a single band, which has the largest transition momentum for the polarization direction of the pump and probe pulses. The chosen band is fitted by a parabola that is shown as the blue dashed line in Fig. 3 (b). The extracted band parameters are as follows: the band gap g is 13.1 eV, and the reduced electron-hole mass µ eh is 1.1m e . Employing the parabolic two-band model with these extracted parameters, we computed the transient optical conductivity based on the pump-probe simulation. As explained above, the parabolic fit is only accurate around the band edge. Therefore, we apply the following attenuation function to the transient optical conductivity to reduce the influence from irrelevant energy regions; Here, g,Γ is the corresponding band gap, which is 13.1 eV. Figure 4 (b) shows the truncated transient conductivity, σ T (ω, T probe )w Γ (ω), with the parabolic two-band model of the Γ point. One can clearly see that the result of the parabolic two-band model in Fig. 4 (b) nicely reproduces the half bottom part of the full ab-initio simulation result in Fig. 2.
To provide the complete description, we combine the results of the two-band models for the M and Γ points by simply summing the transient conductivities in Figs. 4 (a) and (b). The combined transient conductivity is shown in Fig. 5. One sees that the combined result nicely reproduces the full ab-initio result in the whole investigated photon-energy range. Therefore, we can conclude that the fishbone structure in the transient absorption spectrum of monolayer h-BN consists of two major contributions: One is the contribution from the negative-mass electron-hole band at the M point, accounting for the upper half of the fishbone structure. The other is the contribution from the positive-mass band at the Γ point, inducing the lower half of the fishbone structure.
Since the pump-induced interband transitions are completely omitted in the model by construction, the transient absorption of the two-band model in Fig. 5 results only from the IR-induced intraband transitions. Therefore, we clearly demonstrated that the IR-induced intraband transitions play crucial roles and the transient absorption is dominated by the dynamical Franz-Keldysh effect even in this 2D material.
Furthermore, based on the decomposition of the transient absorption spectrum in Fig. 5 into the contributions from each band in Figs. 4 (a) and (b), we found that the upper-half part of the transient absorption spectrum of monolayer h-BN originates from the dynamical Franz-Keldysh effect of a specific band around the M point with the anomalous band dispersion, while the lower-half part originates from the conventional dynamical Franz-Keldysh effect with the positive-mass electron-hole band. As intraband motion is often regarded as semi-classical motion of a free particle, the predominant role of intraband motion in the anomalous dispersion has the potential to induce further interesting dynamics in the 2D material. Finally, we report the transient optical property of monolayer h-BN for a different polarization direction. Here, we set the polarization direction of the pump and probe pulses to the y-axis in Fig. 7. Thus, the polarization direction is parallel to a B-N bond. For this direction, the system does not have the spatial inversion symmetry and it may have even order nonlinear responses. To investigate transient optical properties, we perform pump-probe simulation with the following parameters. For the pump pulse, the peak field strength | E pump | is set to 8.7 × 10 7 V/m, the mean frequency ω pump is set to 1.55 eV/h, and the full pulse-duration T pump is set to 20 fs. For the probe pulse, | E probe | is set to 8.7 × 10 7 V/m, ω probe is set to 15 eV/h, and T probe is set to 1 fs. Under these conditions, we repeat the TDDFT pump-probe simulations by changing the pump-probe time-delay, T delay . Figure 6 (a) shows the computed transient conductivity with the TDDFT pump-probe simulation, and Fig. 6 (b) shows the time profile of the applied pump electric field. In contrast to the 2ω pump -oscillation in Fig. 2, which was the consequence of the inversion symmetry, the transient conductivity shows an oscillatory feature with the frequency of ω pump in time domain. The emergence of the ω pump -oscillation instead of the 2ω pump -oscillation is due to the even-order nonlinear responses, which are forbidden for the x-direction due to inversion symmetry.
Despite of the qualitatively different time-domain structures between Fig. 2 and Fig. 6, both pump-probe configurations show characteristic features in the same photon energy region around 15 eV. According to the above study with the electron-hole pair distribution in Fig. 3 and the two-band model analysis, these features might also be related to the anomalous dispersion of the negative-mass electron-hole band around the M point.
Because the two-band model has spatial inversion symmetry and does not have any even-order nonlinearity, the present model is not appropriate to investigate the ω pump -oscillation in the spectrum of Fig. 6. Therefore, an extension of the model would be required to understand the microscopic mechanism of the ω pump -oscillation structure.  Fig. 1 (a). (b) The time-profile of the applied pump electric field. The result is computed by TDDFT with ALDA.

Summary and conclusions
In this work, we extended the ab-initio pump-probe simulation of ATAS for 2D materials in order to provide a platform to understand the microscopic origin of the features observed in ATAS experiments. As an example, we investigated ultrafast modification of optical properties of monolayer h-BN under intense IR pulses.
First, we performed the pump-probe simulations, setting the polarization of the pump and probe pulses along the x-direction in Fig. 1 (a), which is perpendicular to a B-N bond. As a result of the simulation, we obtained the transient absorption spectrum during the laser irradiation. The resulting spectrum presented an oscillatory behavior with a V-shaped energy dispersion; the fishbone structure. The frequency of the oscillatory behavior was found to be twice the frequency of the pump pulse. In previous studies, similar structures have been observed in bulk materials [9,12,38] and understood in terms of the dynamical Franz-Keldysh effect [39,40].
To clarify the microscopic mechanism of the observed transient absorption spectrum of h-BN, we studied the electronic structure of the material with the DFT and found relevant transitions at the M and Γ points. To elucidate the contribution from these transitions, we further analyzed the transient absorption spectra with a parabolic two-band model. As a result, we identified that the fishbone structure in Fig. 2 consists of two major contributions. One is the dynamical Franz-Keldysh effect around the Γ point. The other is also the dynamical Franz-Keldysh effect but with an anomalous band dispersion around the M point. Therefore, we clearly demonstrated that the IR-induced intraband transitions are significant for the transient absorption even in 2D materials. Furthermore, the discovered mechanism based on intraband motion with the anomalous band dispersion indicates the possibility of introducing exotic dynamical properties of the 2D material since the classical analog of the intraband motion in the anomalous dispersion would be a classical motion of a negative mass particle.
Next, we performed the TDDFT pump-probe simulation, setting the polarization of the pump and probe pulses along the y-direction. For this direction, the inversion symmetry of the material is broken, and the system displays even-order nonlinear responses. Indeed, reflecting the nature of the even-order nonlinearity, the obtained transient absorption spectrum shows an oscillation with the frequency of the pump pulse, ω pump . Therefore, it was demonstrated that the system presents a qualitatively different behavior for different polarization directions. However, despite of the significant differences for the different polarizations in Fig. 2 and Fig. 6, still the characteristic features emerge around the probe photon energy of 15 eV. This fact indicates that the ω pump -oscillation structure in the transient absorption spectrum with the y-direction polarization may also be related to the anomalous band dispersion at the M point. In order to understand the microscopic origin of the ω-oscillation structure, an extension of the parabolic two-band model is required.
The demonstrated ab-initio ATAS simulation for 2D materials is a very general approach that is applicable not only to single-layered materials but also to multi-layered systems including heterostructures. Here, the interplay of multi-layer and interface effects with geometrical tilting, such as moiré patterns, could strongly affect the electron dynamics and be observed in the time domain.
Author Contributions: S.A.S performed the calculations. All the authors participated the interpretation of the results and wrote the paper.