Nonlinear dynamics of electromagnetic field and valley polarization in WSe₂ monolayer

Light–matter interaction is a complex phenomenon that is characterized by two different spatial scales: the micrometer scale for the wavelength of the light and the sub-nanometer scale for the dynamics of electrons.¹ The interaction of laser light with solids is a promising tool to investigate the fundamental physical aspects of material properties in condensed matter physics.^2–4 In recent days, a tremendous amount of research is focused on structural, electronic, and optical properties of thin materials.^5,6 Thin materials of thicknesses less than the wavelength of light display very attractive electronic and optical properties.^7,8 In this regard, two-dimensional (2D) layer systems, such as graphene,^9,10 transitional metal dichalcogenides (TMDs),^11,12 and phosphorene,¹³ are receiving extensive research efforts. Among mentioned 2D materials, TMDs are of particular interest due to broken inversion symmetry and strong spin–orbit coupling (SOC) for valleytronics and spintronics applications. In addition, there is a growing interest to study the nonlinear optical properties of TMDs such as ultrafast carrier dynamics and high order harmonics.^14–17 

The time-dependent density-functional theory (TDDFT) has been proposed as a tool to describe the electron dynamics induced by a time-dependent external potential from the first principles.^18–20 The TDDFT has been most successful in the linear response regime to describe electronic excitation.^21,22 It has also been applied to the nonlinear and nonperturbative dynamics of electrons induced by an intense ultrashort laser pulse.^23–26 The effectiveness of electron dynamics calculations based on the TDDFT is further enhanced by combining it with the Maxwell equations for electromagnetic fields of pulsed light.²⁷ The combined scheme known as the Maxwell-TDDFT formalism has been applied to crystalline solids.^28,29 The Maxwell-TDDFT scheme is quite comprehensive, but a high computational cost limits its applicability. For sufficiently thin materials such as monatomic layers in which macroscopic electromagnetic fields may be treated as spatially uniform within the material, we have developed an efficient approximate method that we call a 2D Maxwell-TDDFT scheme.^30,31 Furthermore, SOC with noncollinear local spin density is incorporated to extend the usefulness of the method to study spin-dependent properties.

In this Letter, we study the linear and nonlinear optical response of the WSe₂ monolayer by using a circularly polarized driving field. The WSe₂ monolayer is chosen to study because of strong SOC and its most promising phenomena of valleytronics. The laser intensity dependence of the valley polarization and light propagation in terms of the transmitted and reflected high harmonics is also investigated by the 2D Maxwell-TDDFT scheme.

Time-dependent calculations are done using an open-source TDDFT package scalable Ab initio light-matter simulator for optics and nanoscience (SALMON). Full details of the SALMON code and its implementation are described elsewhere.^32,33 Here, we briefly describe the 2D Maxwell-TDDFT method for light propagation in thin layers at normal incidence.

We assume that a thin layer is in the xy plane and a light pulse propagates along the z axis. By using the Maxwell equations, we can describe the propagation of the macroscopic electromagnetic fields in the form of the vector potential $A(z,t)$ as

where $J(z,t)$ is the macroscopic current density of the thin layer.

For a very thin layer, we may assume the macroscopic electric field inside the thin layer is spatially uniform. We also approximate the macroscopic electric current density in Eq. (1) as

where $J2D(t)$ is the 2D current density (current per unit area) of the thin layer. We deal with it as a boundary value problem where reflected (transmitted) fields can be determined by the connection conditions at z = 0. From Eq. (2), we obtain the continuity equation of $A(z,t)$ at z = 0 as follows:

where the $A(i), A(r)$⁠, and $A(t)$ are the incident, reflected, and transmitted fields, respectively. From Eqs. (1) and (2), we get the basic equation of the 2D approximation method

Here, $J2D[A(t)](t)$ is the 2D current density that is determined by the vector potential at z = 0, and it is equal to $A(t)(t)$⁠. Time evolution of electron orbitals in a unit cell of the 2D layer driven by $A(t)(t)$ is determined by the time-dependent Kohn–Sham (TDKS) equation. By using the velocity gauge, the TDKS equation for the Bloch orbital $ub,k(r,t)$ (which is a two-component spinor; b is the band index, and $k$ is the 2D crystal momentum of the thin layer) is described as

where the scalar potential $φ(r,t)$ includes the Hartree potential from the electrons and the local part of the ionic pseudopotentials, and we have defined $v̂NLk≡e−ik·rv̂NLeik·r$⁠. Here, $v̂NL$ and $vxc(r,t)$ are the nonlocal part of the ionic pseudopotentials and exchange–correlation potential, respectively. The SOC is incorporated through the j-dependent nonlocal potential $v̂NL$⁠.³⁴ The Bloch orbitals $ub,k(r,t)$ are defined in a box containing the unit cell of the 2D thin layer sandwiched between vacuum regions. The 2D current density $J2D[A(t)](t)$ in Eq. (4) is derived from the Bloch orbitals as follows:

where Ω is the area of the 2D unit cell and the sum is taken over the occupied orbitals in the ground state. In the 2D Maxwell-TDDFT method, coupled Eqs. (4) and (5) are simultaneously solved in real time.

In the weak field limit, we can consider a linear response formalism for the 2D approximation method. The constitutive relation in linear response in the form of the 2D current density $J2D[A(t)](t)$ can be described as follows:

where α and β are the spatial indices and we have introduced $σ̃αβ(t−t′)$ as the 2D electric conductivity of the thin layer. The frequency-dependent 2D conductivity can be obtained by taking the Fourier transformation of $σ̃αβ(t)$⁠. For a weak field limit, we can simply get the relation among the incident, reflection, and transmission by $σ̃αβ(ω)$⁠. For example, for a given incident field $A(i)(t)$⁠, the transmitted field $A(t)(t)$ can be constructed by solving the following equation in frequency space, instead of solving Eq. (4) in time

The WSe₂ monolayer is used to validate our method, and Fig. 1(a) illustrates the crystal structure. The lattice constant is set to $a=b=3.32$ Å. We solve the TDKS equation in the slab geometry with the distance between layers of 20 Å. The dynamics of the 24 (12 electrons for W and 12 electrons for Se) valence electrons are treated explicitly while the effects of the core electrons are considered through norm-conserving pseudopotentials from the OpenMX library.³⁵ The adiabatic local density approximation with Perdew–Zunger functional³⁶ is used for the exchange–correlation. We adopt a spin noncollinear treatment for the exchange–correlation potential.^37,38 The spatial grid sizes and k-points are optimized according to the convergent results. The determined parameter of the grid size is 0.21 Å while the optimized k-mesh is 15 × 15 in the 2D Brillouin zone. The time step size is set to 5 × 10⁻⁴ femtosecond (fs). We consider the circularly polarized laser field with a time profile of $cos4$ envelope shape for the vector potential, which gives the electric field of the applied laser pulse through the following equation:

where $ω$ is the average frequency, $Emax$ is the maximum amplitude of the electric field, and T_P is the pulse duration. We use the frequency of 1.55 eV, T_p is set to 30 fs, and the computation is terminated at 50 fs.

The 2D WSe₂ monolayer has strong SOC, and all the calculations are performed by taking the SOC effect into consideration. The calculated bandgap is ∼1.25 eV. Zeeman-type spin splitting due to SOC is also checked. The value of spin splitting for valence band maxima (VBM) is ∼450 meV, consistent with previous works.^39–41 Figure 1(b) depicts the real and imaginary parts of the optical conductivity $σ̃(ω)$ as a function of the photon energy. $σ̃(ω)$ provides an understanding of the band structure properties, and the linear dependence of $σ̃(ω)$ at small energies with zero value at ω = 0 reveals the pure semiconducting nature of the WSe₂ monolayer. The peaks at higher energies in the optical conductivity correspond to the interband transition from the valence band to the conduction band.

Figure 2 shows a typical electron dynamics calculation of the WSe₂ monolayer. The time profile of the incident circular polarized electric field and the induced electric current at weak and strong intensities is shown in Fig. 2(a). The current at weak intensities is multiplied by a factor of 10 so that the difference between two currents indicates the nonlinear effect at the strong intensity. For a more clear comparison, we have also shown the X–Y view of the current density in Fig. 2(b). At the weak intensity, the temporal evolution of the induced current mostly follows the driving laser profile, indicating that a linear optical response of the WSe₂ monolayer dominates at I = 10¹⁰ W/cm². We also note that there appears substantial current after the incident pulse ends, since the frequency of the applied laser pulse is above the bandgap value. At the strong intensity, the current is initially very close to the case of the weak intensity. However, the current gradually becomes weaker than that expected from the linear response. We consider that this nonlinear effect suppressing the current comes from the saturable absorption that often appears in various 2D materials.^42–44 The saturable absorption takes place by two mechanisms: a decrease in electrons in the valence band by the excitation and an increase in electrons in the conduction band that block the excitation from the valence band. We also note that the current after the incident pulse ends is much weaker than the case of the weak intensity. This indicates that nonlinear effects work to cancel the coherence that produced the delayed current in the linear case.

We show the calculated results of electric excitation energy and the number density of excited electrons in Fig. 3. To see the nonlinear effects, we have scaled up the excitation by a factor of 100 for weak intensity. Note that the weak and strong field curves should coincide and show similar behavior if the response is linear. We first consider the case of the weak intensity. The electronic excitation energy and the number of excited electrons show a similar time profile with a peak around 18 fs and then decrease. The ratio of the two curves is about 1.8 eV and is close to the photon energy of 1.55 eV, confirming the dominance of the one-photon absorption process. The decrease continues even after the incident pulse ends at 30 fs. This decrease is due to emission of light that is described in the present 2D Maxwell-TDDFT scheme and is related to the appearance of the current after the pulse ends in Fig. 2. We next move to the case of the strong intensity. The number of excited electrons is much smaller than estimation by a linear excitation mechanism. This is due to the saturable absorption, as we discussed for the current shown in Fig. 2. However, the electronic excitation energy looks to show mostly linear behavior. We consider that this is due to a cancelation of two effects: nonlinear saturable absorption and multiple absorption of photons. The ratio of the electronic excitation energy to the number of excited electrons is less than 2.0 eV at I = 10¹⁰ W/cm², close to the incident photon energy of 1.55 eV. However, the ratio at I = 10¹² W/cm² is about 3.8 eV that is more than twice the incident photon energy. This indicates the significance of multiple absorption of photons (either by nonlinear multi-photon absorption or multiple absorption of single photons) at I = 10¹² W/cm². The decrease in the number of excited electrons as shown in Fig. 3, despite the multiple absorption of photons is due to the saturable absorption.

To get an understanding of valley pseudospin, we investigate the k-resolved excited electron populations. It is well known that the responses of the conductance for spin-dependent electrons to left- and right-handed resonant circularly polarized light are just the opposite for TMDs.^45–47 Figure 4 shows the excited carrier population by right (σ⁺) and left (σ⁻) circular helicities. A strong circular dichroic effect is evident in the excited state population both at weak and strong intensities. The valley degeneracy is lifted, and pumping with σ⁺ polarized light excites electrons in the $K′$ valley while pumping with σ⁻ polarized light excites the K valley, so-called valley-spin locking. Due to the nonlinear interaction at I = 10¹² W/cm², the excited electrons are not localized around K/$K′$ and start to spread in the Brillouin zone, and an additional peak between the K (K′) and $Γ$-point appears with the nonlinear process. This additional peak induced by the strong field may explain the behavior of the electric current density at I = 10¹² W/cm² as shown in Fig. 2. The excited electron at different points of the Brillouin zone indicates the electronic current with different phases and frequencies. This phase difference decreases the total current after the pulse ends. Overall, the k-resolved excited electron populations demonstrate that the valley polarization state is quite robust against the field strength.

High harmonic generation (HHG) is a promising tool to study strong-field effects and ultrafast electron dynamics in 2D materials. Experimentally, one can primarily measure the fraction of reflected or transmitted pulses of the light, and the 2D Maxwell-TDDFT scheme used for our calculations can describe the propagation of the light wave dynamics as well. Figure 5 shows the HHG spectra including the reflected (⁠$A(r)(t)$⁠) and transmitted (⁠$A(t)(t)$⁠) pulses. Figure 5(a) shows the HHG spectra for a peak intensity of I = 10¹⁰ W/cm² while Fig. 5(b) is for I = 10¹² W/cm². First of all, HHG spectra of the reflected and transmitted pulses completely coincide with each other except for the fundamental frequency, as understood from the continuity equation of Eq. (3). We find that a resonant circularly polarized driver (ω  = 1.55 eV) generates strong high harmonics. For the weak intensity of 10¹⁰ W/cm², the spectrum includes very weak HHG components except for the second harmonic. It was as expected because the electronic current density shows the linear optical response at this intensity. On the other hand, the strong intensity of 10¹² W/cm² [Fig. 5(b)] shows the harmonic order up to 11th order. The WSe₂ monolayer displays both odd and even harmonic peaks. Our results show that the in-plane harmonics of the WSe₂ monolayer under a resonant circular driver has an order of 3n ± 1 (where n is an integer), while the multiple of three harmonic orders are suppressed. The threefold rotational symmetry of the WSe₂ monolayer leads to the selection rule of 3n ± 1. The third harmonic appears in Fig. 5(b) is the exception to this rule, and it may be related to the relative phase effect as explained in Ref. 48.

In this work, we have studied the linear and nonlinear optical response of the WSe₂ monolayer by the classical Maxwell equations combined with TDKS equations to describe the propagation of electromagnetic fields in thin 2D layers. By applying the chiral resonant pulses, the electron dynamics and HHG are shown at weak and strong laser fields. The WSe₂ monolayer shows the linear optical response at the weak intensity. By changing the chirality of the resonant light, a strong circular dichroic effect is observed in the excited state population. A relatively weak laser field shows effective valley polarization while the strong field induces an additional peak of the spin-polarized carrier by the nonlinear process. The WSe₂ monolayer displays the nonlinear optical characteristic of saturable absorption at the strong laser field. The 2D Maxwell-TDDFT scheme offers a clear advantage to measure the HHG in terms of the transmitted and reflected waves. HHG spectra of the reflected and transmitted pulses completely coincide with each other. HHG in the WSe₂ monolayer is observed up to the 11th order, and the circularly polarized resonant driver has a high harmonics of the 3n ± 1 order in the WSe₂ monolayer.

This research was supported by JST-CREST under Grant No. JP-MJCR16N5. This research is also partially supported by JSPS KAKENHI via Grant No. 20H02649 and the MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) under Grant Nos. JPMXS0118068681 and JPMXS0118067246. The numerical calculations are carried out using the computer facilities of the Fugaku through the HPCI System Research Project (Project ID. hp210137), SGI8600 at Japan Atomic Energy Agency (JAEA) and the Multidisciplinary Cooperative Research Program in CCS, University of Tsukuba.