Frequency-resolved Microscopic Current Density Analysis of Linear and Nonlinear Optical Phenomena in Solids

Firstly, we briefly describe the first-principles electron dynamics calculations for solids based on TDDFT Runge and Gross (1984). The details of the method can be referred to fromprevious studies Sato (2021). In the TDDFT based on the Kohn–Sham scheme, the electron dynamics in solids are described by the following time-dependent Kohn–Sham equation:

where $b$ is the band index, and $\boldsymbol{k}$ is the Bloch wavevector. Here, $u_{b\boldsymbol{k}}(\boldsymbol{r},t)$ is a periodic part of the time-dependent Bloch orbital and satisfies the periodic boundary condition, $u_{b\boldsymbol{k}}(\boldsymbol{r},t)=u_{b\boldsymbol{k}}(\boldsymbol{r}+\boldsymbol{a}_{n},t)$, with the lattice vectors, $\boldsymbol{a}_{n}$. The time-dependent Kohn–Sham Hamiltonian, $h_{\boldsymbol{k}}(t)$, is defined as

where $\boldsymbol{A}(t)$ is a vector potential, which is related to a spatially uniform time-dependent electric field as $\boldsymbol{E}(t)=-\frac{1}{c}\dot{\boldsymbol{A}}(t)$. In the present study, the electron–ion interaction $\hat{v}_{\mathrm{ion}}$ is described by the norm-conserving pseudopotential approximation Kleinman and Bylander (1982); Troullier and Martins (1991). Here, $v_{\mathrm{Hxc}}(\boldsymbol{r},t)$ represents the time-dependent Hartree-exchange-correlation potential, and we employ adiabatic local density approximation Perdew and Zunger (1981).

Note that one may add the exchange-correlation vector potential $\boldsymbol{A}_{xc}$ to the vector potential $\boldsymbol{A}(t)$ based on the time-dependent current density functional theory Vignale and Kohn (1996); Vignale et al. (1997); Kootstra et al. (2000). In this work, we set $\boldsymbol{A}_{xc}$ to zero for simplicity.

To describe electron dynamics under specific fields $\boldsymbol{A}(t)$, we solve the time-dependent Kohn–Sham equation, Eq. (1), in the time domain. Note that the initial conditions for the Kohn–Sham orbitals are obtained by solving the static Kohn–Sham equation self-consistently.

where $\epsilon_{b\boldsymbol{k}}$ is the corresponding Kohn-Sham eigenvalue.

Employing the computed time-dependent Kohn–Sham orbitals $u_{b\boldsymbol{k}}(t)$, several physical observables can be evaluated. Among various observables, the macroscopic current density $\boldsymbol{J}(t)$ is an important quantity for studying light–matter interactions, and it is defined as

where $\Omega_{\mathrm{cell}}$ is the unit-cell volume, and $f_{b\boldsymbol{k}}$ is the occupation factor that satisfies

with the number of electrons in te unit-cell, $N_{\mathrm{elec}}$. In Eq. (4), the velocity operator, $\boldsymbol{v}_{\boldsymbol{k}}(t)$, is defined as

where $\boldsymbol{r}$ is the position operator. The velocity operator, $\boldsymbol{v}_{\boldsymbol{k}}(t)$, can be decomposed into a (semi)local part and nonlocal part:

Here, the local and nonlocal parts are defined as follows

Note that the nonlocal contribution originates from the nonlocal part of the pseudopotential approximation, and it vanishes in the all-electron calculations.

Various optical properties of matter can be investigated by analyzing the macroscopic current density $\boldsymbol{J}(t)$ induced by an external electric field $\boldsymbol{E}(t)$. For instance, linear optical properties of solids can be investigated by analyzing the macroscopic current, $\boldsymbol{J}(t)$, induced by a weak perturbation Bertsch et al. (2000). In this work, we evaluate the optical conductivity $\sigma(\omega)$ by applying the Fourier transform to the macroscopic current density $\boldsymbol{J}$ that is induced by an impulsive distortion, $\boldsymbol{E}(t)=\boldsymbol{e}_{z}E_{0}\delta(t)$, as

where $\gamma$ is a phenomenological damping parameter, and we set $\gamma$ to $0.1$ eV.

Similarly, nonlinear optical properties can be investigated from the current $\boldsymbol{J}(t)$ induced by an intense laser field Goncharov (2013); Grüning et al. (2016); Uemoto et al. (2019). Furthermore, by employing the pump–probe setup, one can investigate transient optical properties of solids in the presence of laser fields in the time domain Sato et al. (2014). Hence, the macroscopic current density $\boldsymbol{J}(t)$ plays a crucial role in the investigation of the optical properties of solids.

We perform microscopic current density analysis in the frequency domain to gain microscopic insight into light-induced phenomena. Here, we consider the following microscopic current density

Note that the local current density cannot be straightforwardly defined for the nonlocal part of the velocity operator, $\boldsymbol{v}^{\mathrm{nonlocal}}_{\boldsymbol{k}}(t)$. Furthermore, the contribution from the nonlocal operator to the macroscopic current density, $\boldsymbol{J}(t)$, vanishes in the all-electron calculations. Hence, in this study, we consider the contribution only from the local part of the velocity operator, $\boldsymbol{v}^{\mathrm{local}}_{\boldsymbol{k}}(t)$, to define the microscopic current density, $\boldsymbol{j}(\boldsymbol{r},t)$, in Eq. (11). When the contribution of the nonlocal operator, $\boldsymbol{v}^{\mathrm{nonlocal}}_{\boldsymbol{k}}(t)$, is negligible, the spatial average of the microscopic current density, $\boldsymbol{j}(\boldsymbol{r},t)$, in the unit-cell reproduces the macroscopic current density.

Hence, the microscopic current density, $\boldsymbol{j}(\boldsymbol{r},t)$, is expected to provide microscopic insight into light-induced phenomena in solids.

Having established the connection between the optical properties of solids and microscopic current density via the macroscopic current density, we further analyze the microscopic current density in the frequency domain. For this analysis, we introduce frequency-resolved microscopic current dneisty, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega)$:

where $W(t)$ is a window function for the reduction of the numerical error due to the finite simulation time. For practical calculations in this study, we employ the following window function:

in the duration $-T_{\mathrm{pulse}}/2<t<T_{\mathrm{pulse}}/2$ and zero outside. Here, $T_{\mathrm{pulse}}$ is the full duration of applied laser pulses, which will be introduced later along with the pulse shape of the applied fields.

By analyzing the frequency-resolved microscopic current density, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega)$, in Eq. (13), we can gain atomic-scale insight into linear and nonlinear optical phenomena in solids at a specific frequency of optical responses. For instance, we can study whether localized electrons on specific bonds or atoms are responsible for a certain optical response or whether delocalized electrons play a central role by visualizing the frequency-resolved microscopic current density, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega)$.

The discrepancy between the Fourier transform of current density, $\boldsymbol{j}(\boldsymbol{r},t)$, and that of electron density, $\rho(\boldsymbol{r},t)$, from the perspective of optical responses of solids is noteworthy. If the nonlocal part of the velocity operator is negligible, the microscopic current density and the electron density are related via the continuity equation.

Although certain information pertaining to the microscopic current density, $\boldsymbol{j}(\boldsymbol{r},t)$, can be reproduced based on the electron density, $\rho(\boldsymbol{r},t)$, via Eq. (15), the divergence-free part of $\boldsymbol{j}(\boldsymbol{r},t)$ cannot be reproduced, since such component does not contribute to the continuity equation. Therefore, when the divergence-free component of $\boldsymbol{j}(\boldsymbol{r},t)$ plays a crucial role in optical response, the electron density may not contain sufficient information to describe such response. For instance, let us consider an ideal noninteracting Fermi gas as an approximated electronic system. The ground state electron density of such a system is constant in real-space and constitutes a sphere in momentum space, the so-called Fermi sphere. We consider applying a homogeneous electric field to such a system. Owing to the field application, the system exhibits homogeneous current density in realsspace due to the shift of the Fermi sphere in momentum space, while the electron density remains constant in both space and time. Hence, the electron density does not reflect the light-induced electron dynamics of ideal Fermi gases. This example indicates that electron density analysis is not sufficient for analyzing the optical responses of delocalized electrons, such as free electrons in metals, although density analysis is useful for isolated systems. In contrast, the microscopic current density can be used for both localized and delocalized electrons because it is directly linked to the macroscopic current density $\boldsymbol{J}(t)$ and optical responses via the Maxwell equation. Hence, the microscopic current density affords a more comprehensive description of the optical responses of matter.

First, we examine the microscopic current density analysis based on linear responses of a simple metal, aluminum. For the description of bulk aluminum, we employ a cubic unit cell that contains four aluminum atoms. The aluminum unit cell is discretized into $16^{3}$ real-space grid points, and the corresponding first Brillouin is also discretized into $24^{3}$ $k$-points. Aluminum atoms are described by a norm-conserving pseudopotential method, treating $3s$ and $3p$ electrons as valence electrons.

To gain insight into electron dynamics in simple metals, we revisit the electron density in aluminum. Figure 1 shows the computed electron density, $\rho(\boldsymbol{r})$, in the ground state of aluminum. Electrons are delocalized in aluminum. These delocalized electrons correspond to conduction electrons, and they are expected to play a central role in the optical responses of aluminum.

Figure 2 shows the optical conductivity $\sigma(\omega)$ of bulk aluminum computed via TDDFT calculations. Note that the conductivity and dielectric function, $\epsilon(\omega)$, are connected via $\epsilon(\omega)=1+\frac{4\pi i}{\omega}\sigma(\omega)$, and the photoabsorption is closely related to the real part of the conductivity, $\mathrm{Re}\left[\sigma(\omega)\right]$, and the imaginary part of the dielectric function, $\mathrm{Im}\left[\epsilon(\omega)\right]$. As shown in Fig. 2, both real and imaginary parts of the conductivity show positive values in the investigated photon energy range. This is a typical metallic response with respect to optical fields below the plasma frequency, and the Drude model can describe it suitably.

To investigate the metallic response of aluminum, we compute electron dynamics under a laser pulse described by the following vector potential,

in the duration $-T_{\mathrm{pulse}}/2<t<T_{\mathrm{pulse}}/2$ and zero outside. Here, $\boldsymbol{e}_{z}$ is the unit vector along $z$-axis, which is $(001)$-direction of the unit cell, $E_{0}$ is the peak field strength, $T_{\mathrm{pulse}}$ is the full pulse duration, and $\omega_{L}$ is the mean frequency of the laser field. For the calculations of electron dynamics in aluminum, we set $T_{\mathrm{pulse}}$ to $20$ fs, $\omega_{L}$ to $1.55$ ev/$\hbar$, and $E_{0}$ to $2.75\times 10^{8}$ V/m. Employing the computed time-dependent current density, $\boldsymbol{j}(\boldsymbol{r},t)$, we further evaluate the frequency-resolved microscopic current density using Eq. (13) based on the electron dynamics calculation.

Figure 3 shows the $z$-component of the frequency-resolved microscopic current density, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega=\omega_{L})$, at the driving laser frequency $\omega_{L}$. Since the carrier wave of the electric field $\boldsymbol{E}(t)$ is derived from the vector potential, Eq (16) and it is proportional to $\cos(\omega_{L}t)$, the real part of $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega_{L})$ corresponds to the real part of the conductivity in Fig. 2, while the imaginary part of $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega_{L})$ corresponds to the imaginary part of the conductivity.

As observed in Fig. (3), both real and imaginary parts of $\tilde{j}_{z}(\boldsymbol{r},\omega_{L})$ in aluminum exhibit quasi-homogeneous positive contributions to the whole space. This indicates that the optical responses of simple metals are dominated by delocalized free electrons in the conduction bands, and the frequency-resolved microscopic current density suitably captures this metallic nature of the optical responses with free carriers.

We subsequently perform Fourier analysis on the microscopic current density with linear responses of a typical semiconductor, i.e., silicon. For the description of silicon, we employ a cubic unit cell that contains eight silicon atoms. The silicon unit cell is discretized into $20^{3}$ real-space grid points, and the corresponding first Brillouin zone is also discretized into $24^{3}$ $k$-points. Silicon atoms are described via a norm-conserving pseudopotential method, treating $3s$ and $3p$ electrons as valence electrons.

To obtain microscopic insight into the electronic properties of silicon, we revisit the electron density of silicon. Figure 4 shows the computed electron density $\rho(\boldsymbol{r})$ of silicon in the ground state. The electron density is concentrated between the two silicon atoms, forming covalent bonds. In contrast to delocalized electrons in aluminum in Fig. 1, electrons are highly localized in silicon. Hence, dielectric responses of silicon are expected to be dominated by the dynamics of localized electrons at covalent bonds.

Here, we elucidate the dielectric responses of silicon. For this purpose, we compute the optical conductivity of silicon. Figure 5 shows the real and imaginary parts of the optical conductivity, $\sigma(\omega)$, of silicon computed through TDDFT. The computed optical gap of silicon is approximately $2.4$ eV, and the real part of the conductivity almost vanishes below the optical gap, reflecting the absence of photoabsorption below the gap. In contrast to the aluminum case in Fig. 2, the imaginary part of the conductivity of silicon becomes negative below the optical gap. The positive value of the imaginary part of the conductivity of aluminum in the low-frequency region, i.e., at $1.55$ eV, can be explained in terms of the response of free carriers using the Drude model, whereas the negative value of the conductivity of silicon in the low-frequency regime can be explained in terms of the response of bound electrons using the Lorentz model, which is a classical model based on harmonic oscillators. This also indicates that the optical responses of silicon in the low photon energy region below the optical gap are dominated by the bound electrons at the covalent bonds.

Note that the computed optical gap of silicon, $2.4$ eV, underestimates the experimental value of $3.4$ eV Aspnes and Studna (1983). This underestimation is caused by the local density approximation used in the present study. Nevertheless, this underestimation of the gap does not affect the conclusions of this research because the scope of this study is to establish a scheme for analyzing the optical responses of solids and demonstrate its feasibility. Note that the proposed microscopic current density analysis can be combined with any theoretical methods as long as the microscopic current density $\boldsymbol{j}(\boldsymbol{r},t)$ is provided. Hence, the frequency-resolved current density analysis can be easily extended to TDDFT calculations with more advanced exchange-correlation functionals, which correct the gap underestimation problems. Furthermore, one may apply the microscopic current density analysis even to the time-dependent current density functional theory and the exact many-body Schrödinger equation.

To develop microscopic insight into the optical response of silicon, we compute the electron dynamics under the field described with Eq. (16). Here, we use the same parameters as the analysis of aluminum in Sec. III.1. We subsequently evaluate the frequency-resolved microscopic electron density at the frequency of the laser field as $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega=\omega_{L})$ with $\omega_{L}=1.55$ eV/$\hbar$. Figure 6 shows the $z$-component of $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega_{L})$. Consistent with the optical conductivity at $\omega_{L}=1.55$ eV/$\hbar$, the real part of the frequency-resolved microscopic current density is almost zero in Fig 6 (a), while the imaginary part shows negative values. Furthermore, by comparing Fig. 6 (b) with Fig. 4, the microscopic current density of silicon is found to be concentrated in the covalent bond region. Hence, the frequency-resolved microscopic current density analysis indicates that the optical responses of silicon below the optical gap are dominated by the responses of the bound electrons localized at the covalent bonds.

Having demonstrated that the frequency-resolved microscopic current density analysis can capture both metallic responses of delocalized free carriers and dielectric responses of localized electrons at covalent bonds, we subsequently apply Fourier analysis to nonlinear optical responses. As an example for the demonstration, we consider the third harmonic generation from diamond. For the description of diamond, we employ a cubic unit cell that contains eight carbon atoms. The diamond unit cell is discretized into $20^{3}$ real-space grid points, and the corresponding first Brillouin zone is also discretized into $24^{3}$ $k$-points. Carbon atoms are described by a norm-conserving pseudopotential method, treating $2s$ and $2p$ electrons as valence electrons.

To gain microscopic insight into the electronic properties of diamond, we revisit the electron density. Figure 7 shows the ground state electron density in diamond. The covalent bonds are found to be formed between carbon atoms. Although both diamond and silicon show clear covalent bonds, electrons are more localized around the atoms than the bond regions in the case of diamond, whereas more significant electron localization can be found around the bond regions in the case of silicon. Hence, compared with the case of silicon, localized electron dynamics around carbon atoms are expected to be more important in diamond rather than the bond region.

Before investigating the nonlinear optical responses of diamond, we revisit the linear optical properties. Figure 8 shows the optical conductivity of diamond computed via TDDFT. The computed optical gap of diamond is approximately $5.5$ eV, which is underestimated from the experimental gap of $7.1$ eV Logothetidis et al. (1992), and the real part of the conductivity is almost zero below the optical gap. In contrast to the real part, the imaginary part shows finite negative values below the optical gap. These behaviors are consistent with those of silicon in Fig. 5.

As an example of the nonlinear optical effects, here we investigate the third harmonic generation from diamond. For this purpose, we compute the electron dynamics under a vector potential described by Eq. (16) by setting $\omega_{L}$ to $1.0$ eV/$\hbar$, $E_{0}$ to $8.7\times 10^{8}$ V/m, and $T_{\mathrm{pulse}}$ to $20$ fs. Note that the corresponding three-photon energy is $3\omega_{L}\hbar=3.0$ eV, and it is below the optical gap of $5.5$ eV. Hence, the three-phonon process is still an off-resonant process under the present condition.

Figure 9 shows the $z$-component of the frequency-resolved microscopic current density, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega=3\omega_{L})$, in diamond at the frequency of $3\omega_{L}$, which is the three times of the driving laser frequency. Hence, the microscopic current density in Fig. 9 reflects the nature of the third-order harmonic generation. Similar to the cases of the linear responses of aluminum in Sec. III.1 and silicon in Sec. III.2, the real and imaginary parts of the frequency-resolved microscopic current density, $\tilde{\boldsymbol{j}}(\boldsymbol{r},\omega)$, at $\omega=3\omega_{L}$ correspond to the real and imaginary parts of the nonlinear conductivity tensors, $\sigma^{(3)}_{ijkl}$, respectively. Reflecting the fact that the real part of the nonlinear conductivity is zero due to the off-resonant condition, the real part of the frequency-resolved microscopic current density in Fig. 9 is almost zero in the entire space. In contrast, the imaginary part shows the localized structures around diamond atoms. This indicates that the third harmonic generation in diamond in the off-resonant regime is caused by the dynamics of bound electrons around carbon atoms rather than the covalent bond region. Furthermore, more microscopic insight is obtained: a positive contribution (red) appears around the relatively outer region of the bound electron dynamics, while a negative contribution (blue) appears around the inner region. Hence, there is an internal cancellation in the third-order harmonic generation in diamonds. Recently, a similar cancellation but in $k$-space has been discussed for the high-order harmonic generation in graphene Sato et al. (2021). Furthermore, the enhancement of the high-order harmonic generation has been proposed by suppressing the cancellation with an external degree of freedom. By extending the knowledge of the previous work, one may further enhance the third-order harmonic generation by suppressing the cancellation with an external degree of freedom, such as the secondary light field.

At the end of this section, we discuss a contribution from the nonlocal part of pseudopotentials to the macroscopic current density. The microscopic current density in Eq. (11) is defined only by the local part of the velocity operator in Eq. (8), ignoring the nonlocal part in Eq. (9). Hence, the spatial average of the microscopic current may not fully reproduce the macroscopic current density $\boldsymbol{J}(t)$ if the nonlocal part of the velocity operator exists due to the pseudopotential approximation. Nevertheless, the contribution from the nonlocal part of pseudopotentials does not prevent the extraction of microscopic insight at the atomic scale. Because the nonlocal part of pseudopotentials is constructed with atomic orbitals and is localized within a cutoff radius, the contribution from the nonlocal part to the macroscopic current density is highly localized at atoms, indicating that the nonlocal contribution can be assigned to highly localized atomic responses.

Practically, the contribution from the nonlocal part of the velocity operator to the macroscopic current is often small. To assess the significance of the nonlocal contribution, we compute the power spectrum of emitted light from diamond by applying the Fourier transform to the macroscopic current density obtained by the above calculation as

Figure 10 shows the computed power spectrum, $I(\omega)$. The red solid line shows the result obtained using the macroscopic current density with the nonlocal part of the velocity operator, whereas the blue dashed line shows the result without the nonlocal contribution. As seen from Fig. 10, both responses at the fundamental frequency $\omega=1$ eV/$\hbar$ and the third harmonic frequency $\omega=3$ eV/$\hbar$ are well described, even without the contribution from the nonlocal velocity operator.

We evaluate the following quantity to quantify the error due to the lack of nonlocal contribution:

where $\boldsymbol{J}_{\mathrm{w/~{}nonloc}}(t)$ is the macroscopic current density computed by including the contribution from the nonlocal part of the velocity operator, and $\boldsymbol{J}_{\mathrm{w/o~{}nonloc}}(t)$ is the macroscopic current density computed by excluding the nonlocal contribution. The computed error $I_{\mathrm{Error}}$ is $8.1\times 10^{-2}$ at $\omega=3\omega_{L}$. Hence, in the present analysis of the third harmonic generation, we can safely ignore the nonlocal contribution in the velocity operator.