Momentum-resolved TDDFT algorithm in atomic basis for real time tracking of electronic excitation

TL;DR

This paper introduces a momentum-resolved real-time TDDFT algorithm using atomic basis functions, enabling efficient simulation of ultrafast electronic excitations in extended materials like graphene.

Contribution

The authors develop a novel momentum-resolved rt-TDDFT method with atomic basis and implement electromagnetic field gauges, advancing the simulation of ultrafast dynamics in solids.

Findings

Successfully simulated elementary excitations in graphene

Observed excitation modes distinguishable in momentum space

Demonstrated computational efficiency for extended systems

Abstract

Ultrafast electronic dynamics in solids lies at the core of modern condensed matter and materials physics. To build up a practical ab initio method for studying solids under photoexcitation, we develop a momentum-resolved real-time time dependent density functional theory (rt-TDDFT)algorithm using numerical atomic basis, together with the implementation of both the length and vector gauge of the electromagnetic field. When applied to simulate elementary excitations in two-dimensional materials such as graphene, different excitation modes, only distinguishable in momentum space, are observed. The momentum-resolved rt-TDDFT is important and computationally efficient for the study of ultrafast dynamics in extended systems.

Figures

Equations

$$S_{i\alpha ,j\beta ,\mathbf{k}}=\sum _s{e}^{-i\mathbf{k}\cdot {\mathbf{R}}_s}⟨{\xi }_{i\alpha }(\mathbf{r}+{\mathbf{R}}_s+{\mathbf{b}}_i)⟩{\xi }_{j\beta }(\mathbf{r}+{\mathbf{b}}_j),$$

$$H_{i\alpha ,j\beta ,\mathbf{k}}=\sum _s{e}^{-i\mathbf{k}\cdot {\mathbf{R}}_s}⟨{\xi }_{i\alpha }(\mathbf{r}+{\mathbf{R}}_s+{\mathbf{b}}_i)⟩\hat{H}{\xi }_{j\beta }(\mathbf{r}+{\mathbf{b}}_j),$$

$$\begin{array}{rl}\hat{H}=& \hat{T}+\sum V_I^{local}(\mathbf{r})+\sum V_I^{KB}+V^H(\mathbf{r},\rho (\mathbf{r}))\\ & +V^{XC}(\mathbf{r},\rho (\mathbf{r}))+V^{ext}(\mathbf{r})\end{array}$$

$$V^{ext}(\mathbf{r},t)=-\mathbf{E}(t)\cdot \mathbf{r}.$$

$$\mathbf{E}(t)={\mathbf{E}}_0\cos \left(2\pi ft+\varphi \right)\exp \left[-\frac{(t-t_0{)}^2}{2{\sigma }^2}\right],$$

$$V^{ext}(\mathbf{r}+{\mathbf{R}}_s,t)=-\mathbf{E}(t)\cdot (\mathbf{r}+{\mathbf{R}}_s)\ne -\mathbf{E}(t)\cdot \mathbf{r}.$$

$$E_{\mu }(\mathbf{r},t)=\{\begin{array}{ll}{E}_{\mu }(t)& ϵ<x_{\mu }<L_{\mu }-ϵ,\\ -E_{\mu }(t)L_{\mu }/2ϵ& -ϵ<x_{\mu }<+ϵ.\end{array}$$

$$\mathbf{E}=-\frac{1}{c}\frac{\partial \mathbf{A}}{\partial t};\mathbf{A}=-c\int \mathbf{E}dt.$$

$$H=\frac{1}{2m}(\hslash \mathbf{k}-\frac{e}{c}\mathbf{A}{)}^2=\frac{1}{2m}(\hslash k+e\int \mathbf{E}dt{)}^2=\frac{{\hslash }^2}{2m}(\mathbf{k}+{\mathbf{k}}_{\mathbf{A}}{)}^2,$$

$${\mathbf{k}}_{\mathbf{A}}=\frac{e}{\hslash }\int \mathbf{E}dt=\sqrt{2}\int \mathbf{E}dt.$$

$$|u_{n\mathbf{k}}(\mathbf{r},t_2)⟩=\exp \left[-i{S}_{\mathbf{k}}^{-1}(t^{\prime })H_{\mathbf{k}}(t^{\prime })\Delta t\right]|u_{n\mathbf{k}}(\mathbf{r},t_1)⟩.$$

$$|u_{n\mathbf{k}}(\mathbf{r},t_2)⟩=\frac{1-i{S}_{\mathbf{k}}^{-1}{H}_{\mathbf{k}}\Delta t/2}{1+i{S}_{\mathbf{k}}^{-1}{H}_{\mathbf{k}}\Delta t/2}|u_{n\mathbf{k}}(\mathbf{r},t_1)⟩.$$

$$\begin{array}{rl}{\rho }_{i\alpha ,j\beta }(t_2)& =\sum _n\sum _{\mathbf{k}}{q}_{n,\mathbf{k}}|u_{n\mathbf{k}}(\mathbf{r},t_2)⟩⟨u_{n\mathbf{k}}(\mathbf{r},t_2)|\\ & =\sum _n\sum _{\mathbf{k}}{q}_{n,\mathbf{k}}{c}_{n,i\alpha ,\mathbf{k}}^{\ast }(t_2)c_{n,j\beta ,\mathbf{k}}(t_2),\end{array}$$

$$|u_{n\mathbf{k}}(\mathbf{r},t_2)⟩=\sum _{j\beta }{c}_{n,j\beta ,\mathbf{k}}(t_2){\xi }_{j\beta }(\mathbf{r}).$$

$$\max \left\{\left|{\rho }_{i\alpha ,j\beta }^{new}-{\rho }_{i\alpha ,j\beta }\right|\right\}<\eta ,$$

$$\rho =(1-w)\rho +w{\rho }_{new},$$

$$P_{n{n}^{\prime }\mathbf{k}}={\left|C_{n{n}^{\prime }\mathbf{k}}\right|}^2={\left|⟨v_{n\mathbf{k}}⟩S_k{u}_{n^{\prime }\mathbf{k}}\right|}^2,$$

$$H_{\mathbf{k}}|v_{n\mathbf{k}}(\mathbf{r})⟩=E_{n\mathbf{k}}{S}_{\mathbf{k}}|v_{n\mathbf{k}}(\mathbf{r})⟩.$$

$$q_{n\mathbf{k}}=\sum _{n^{\prime }\in n_{\mathbf{k},occ}}{q}_{n^{\prime }\mathbf{k}}{P}_{n{n}^{\prime }\mathbf{k}},$$

$$\mathbf{j}=-i\frac{e\hslash }{m}\sum _n(⟨u_{n\mathbf{k}}⟩\nabla u_{n\mathbf{k}}-⟨u_{n\mathbf{k}}⟩\nabla {u_{n\mathbf{k}}}^{\ast }),$$

$$\Delta =\frac{1}{T}{\int }_0^T|E_{ex}(t)-E_{ex}^{ref}(t)|dt,$$

$$E_{ex}(t)=E_{KS}(t)-E_{KS}(t=0),$$

$$P_{\mu }(\omega )={\alpha }_{\mu ,\nu }(\omega )E_{\nu }(\omega ).$$

$$\int P_{\mu }(t)\exp (i\omega t)dt={\alpha }_{\mu ,\nu }(\omega )\int E_{\nu }(t)\exp (i\omega t)dt.$$

$${\alpha }_{\mu ,\nu }(\omega )=\frac{\int P_{\mu }(t)\exp (i\omega t)dt}{\int E_{\nu }(t)\exp (i\omega t)dt}.$$

$$E_{\nu }^{\theta }(\omega )=\int E_{\nu 0}[1-\theta (t)]\exp (i\omega t)dt=\frac{E_{\nu }^0}{i\omega }.$$

$$E_{\nu }^{\theta }(t)=E_{\nu }^0[1-\theta (t)]=\{\begin{array}{ll}{E}_{\nu }^0& t\le 0\\ 0& t>0\end{array},$$

$${\alpha }_{\mu ,\nu }(\omega )=\frac{i\omega }{E_{\nu }^0}\int P_{\mu }(t)\exp (i\omega t)dt.$$

$$n(t)=\sum _{unocc}{q}_{n\mathbf{k}}(t),$$

$$q_{n{K}_1}(t)=A_0\cos \left(\frac{{\omega }_{K_2}-{\omega }_r}{2}t\right)\cos \left(\frac{{\omega }_{K_2}+{\omega }_r}{2}t\right).$$

Full text

I Introduction

Real-time (rt) time dependent density functional theory (TDDFT) is an efficient ab initio method to study electron dynamics in complex electron-nuclear systems in both the ground state and excited state. Compared with other widely used approaches such as frequency domain TDDFT, quasi-particle GW, and Bethe-Salpeter equations, rt-TDDFT has two major advantages: (i) Time-dependent Kohn-Sham (TDKS) equations in rt-TDDFT include all nonlinear effects and are intrinsically non-perturbative, making rt-TDDFT a better tool to describe materials in a strong field and (ii) rt-TDDFT directly provides complete information on real time evolution of electronic wavefunctions together with ionic movements, presenting a unique way for real-time tracking ultrafast dynamics and complex phenomena far from equilibrium. Thus, rt-TDDFT is a natural choice for the exploration of strong field physics and ultrafast phenomena. Motivated by the rapid developments in ultrafast experimental techniques, e.g., attosecond based spectroscopy Krausz and Ivanov (2009), ultrastrong laser sources Kruchinin, Krausz, and Yakovlev (2018) and free electron X-ray lasers Pellegrini, Marinelli, and Reiche (2016), rt-TDDFT is drawing more and more attention as a method to simulate ultrafast phenomena in the current line of research frontiers.

Nevertheless, rt-TDDFT is not widely used as the method of choice in the literature, being much less popular than other density functional theory (DFT) based approaches such as $\Delta$SCF, DFT+U, frequency-domain TDDFT, etc. Thus, numerical atomic orbitals (NAO) have been a common choice to dramatically reduce computation cost for simulating complex materials and have been widely used in DFT codes such as SIESTA Soler et al. (2002); Ordejón, Artacho, and Soler (1996) and OpenMX Ozaki (2003) and rt-TDDFT implementations by A. Tsolakidis Tsolakidis, Sánchez-Portal, and Martin (2002) and X. Li Li et al. (2005a); Isborn, Li, and Tully (2007a). The biggest advantage of using NAO is the extremely small computational cost. To describe a system with $N_{a}$ atoms, only about $10\times N_{a}$ NAOs are required, while $10^{3}-10^{4}\times N_{a}$ real space grids or plane waves have to be invoked. In addition, with a relatively small real-space cutoff for NAOs, the order-$N$ linear scaling with respect to system size can be achieved. Since a major difficulty in developing rt-TDDFT is its extreme time consumption due to the use of ultrasmall time step (on the order of $\sim$1 attosecond), NAO based $\mathbf{k}$-resolved rt-TDDFT is very promising for simulating realistic condensed matter systems, complex materials, and interfaces with a long simulation time.

Most previous rt-TDDFT investigations focus on the photoabsorption and related properties of finite-size zero-dimensional (0D) systems (atoms/molecules/nanoparticles) including optical spectra, Yamamoto, Noguchi, and Watanabe (2006); Qian et al. (2006); Tong and Chu (2001, 1998, 1997); Heslar et al. (2007); Nobusada and Yabana (2004); Ganeev et al. (2015); Lopata et al. (2012); Fernando, Balhoff, and Lopata (2015); Tussupbayev et al. (2015); Lopata and Govind (2013); Raghunathan and Nest (2012); Williams-Young, Goings, and Li (2016); Bruner, LaMaster, and Lopata (2016); Provorse, Habenicht, and Isborn (2015); Fischer, Cramer, and Govind (2015); Repisky et al. (2015); Lopata and Govind (2011); Nguyen and Parkhill (2015); Zheng et al. (2016), excited state dynamics Donati et al. (2016); Petrone et al. (2014); Chapman, Liang, and Li (2011), solvation effect Donati et al. (2017a); Ding et al. (2015a); Chapman, Liang, and Li (2013); Nguyen et al. (2012); Liang et al. (2012), relativistic effect variationally Kasper et al. (2018); Goings et al. (2016), photochemical stability Haruyama, Hu, and Watanabe (2012); Haruyama et al. (2012); Hu et al. (2013); Silaeva et al. (2015); Yan, Wang, and Meng (2016), and recently plasmonic excitations Ding et al. (2014); Donati et al. (2017b, 2018); Manjavacas et al. (2014); Barbry et al. (2015); Townsend and Bryant (2012); Ma, Wang, and Wang (2015); Yan, Jacobsen, and Thygesen (2011); Song et al. (2012); Yan, Yuan, and Gao (2007); Gao and Yuan (2011); Gao (2015); Song, Nordlander, and Gao (2011). In 0D systems, only single $\Gamma$ point is needed in the reciprocal space sampling. Thus, the $\Gamma$-only algorithm is overwhelming, as commonly implemented and used in the majority of rt-TDDFT simulations. However, to study photoexcitation and electronic dynamics in extended systems, $\Gamma$-only k-point sampling is insufficient and momentum-resolved ($\mathbf{k}$-resolved) sampling in the reciprocal space is required.

An important advantage of using $\mathbf{k}$-resolved rt-TDDFT is computational efficiency. With ${\Gamma}$-only TDDFT, to get the accurate charge density and ionic forces, an extraordinary large supercell has to be invoked. Many previous studies on extended systems belong to this scenario, Miyamoto, Rubio, and Tománek (2006); Miyamoto (2007a); Krasheninnikov, Miyamoto, and Tománek (2007); Miyamoto (2007b); Zhang and Miyamoto (2009); Miyamoto, Zhang, and Tom??nek (2010); Zhang and Miyamoto (2012); Zhang, Miyamoto, and Rubio (2012); Miyamoto et al. (2013) including our recent studies on ultrafast electron-hole dynamics in dye-sensitized solar cells, Meng and Kaxiras (2010); Meng, Ren, and Kaxiras (2008); Ma, Jiao, and Meng (2014, 2013); Jiao, Ding, and Meng (2011); Jiao, Ma, and Meng (2013) charge separation in van der Waals heterojunctions, Zhang et al. (2017) and nonthermal melting of silicon. Lian, Zhang, and Meng (2016) Using $\mathbf{k}$-resolved algorithms, and at the same accuracy level, the supercell size as well as the computational cost, can be largely reduced, as will be demonstrated later.

Besides technical advantages, $\mathbf{k}$-resolved algorithm introduces the important $\mathbf{k}$-space resolution and a new degree of freedom, which is essential to describe key quantities and important physics in condensed matter materials such as time-dependent band structures, quasiparticles, and valley dynamics. Only rt-TDDFT with $\mathbf{k}$-resolved sampling can provide essential information concerning the real time evolution of material properties.

Although $\mathbf{k}$-resolved rt-TDDFT algorithms have been implemented by several groups Bertsch et al. (2000); Marques (2003); Castro et al. (2006); Andrade et al. (2015) and applied for both semiconductors, Sato et al. (2015a, b); Otobe et al. (2016); Yabana et al. (2012); Shinohara et al. (2010a, b); Otobe, Yabana, and Iwata (2009); Otobe et al. (2008); Wachter et al. (2014) and metals, Krieger et al. (2015); Elliott et al. (2016); Schleife et al. (2012); Yost, Yao, and Kanai (2017) these implementations employ either real space grids or planewaves as basis sets. With a much smaller basis set, the implementation of $\mathbf{k}$-resolved rt-TDDFT algorithms with NAO basis has advantages in efficiency. To take the advantages of NAOs, a new framework and a more complicated implementation of rt-TDDFT are required.

In this work, we strive to tackle the major challenges mentioned above in NAO-based rt-TDDFT. We have successfully developed the $\mathbf{k}$-resolved rt-TDDFT algorithm based on local atomic basis sets using numerical atomic orbitals. Both the length and vector gauge of the electromagnetic field have been implemented. This approach enables rt-TDDFT calculations of solids and surfaces using rather simple unit cells, reducing computational cost by several orders of magnitudes. Moreover, momentum-resolved electron dynamics in the excited states can be tackled by this approach. For instance, $\mathbf{k}$ selective photoexcitations in graphene are demonstrated here, where three distinct photoexcitation modes located at different $\mathbf{k}$points in the reciprocal space are induced upon laser illumination. This kind of $\mathbf{k}$-dependent electronic dynamics is ubiquitous in extended systems such as periodic solids and interfaces. Therefore, we expect highly efficient $\mathbf{k}$-resolved rt-TDDFT algorithms employing local bases be an important development and will be widely used in first-principles simulations of ultrafast phenomena under a strong field and optimal control of quantum materials.

II Methodology

The main framework of $\mathbf{k}$-resolved rt-TDDFT algorithm is inherited from an earlier single-$\Gamma$ version of Time Dependent ab initio Package (TDAP), Meng and Kaxiras (2008) which is based on the SIESTA Soler et al. (2002); Ordejón, Artacho, and Soler (1996) package. In such a rt-TDDFT algorithm, the flowchart of a given ionic step is shown in Fig. 1. Each process is described in detail in the Secs. II A-II G, marked with the same labels as in Fig. 1. Here atomic units $\hbar=m_{e}=e=1$ are used throughout this work.

II.1 Hamiltonian and overlap matrix

Adopting periodical boundary conditions, the lattice of an extended system are denoted as $\mathbf{R}_{s}$ ($s=1,2,3,...,N$) and the atoms $i$ in the unit cell are located at positions $\mathbf{b}_{i}$, where $N$ is truncated to construct a finite supercell. A set of numerical atomic-centered orbitals (NAOs) $\{\xi_{i\alpha}\}$ is associated with each atom in the simulated system, where $\alpha$ denotes both the orbital and angular quantum number of an atomic orbital, each expressed in multiple radial basis functions $\zeta$ Soler et al. (2002). Here, since all the operators and functions are time-dependent, we only denote the explicit dependence on $t$ as $f(t)$ and omit $t$ for implicit dependence.

Overlap matrix $S_{\mathbf{k}}$ and Hamiltonian $H_{\mathbf{k}}$ at the each $\mathbf{k}$ point in the reciprocal space are expressed with NAOs:

[TABLE]

[TABLE]

where

[TABLE]

is the Hamiltonian operator. Here $\hat{T}=\frac{1}{2}\nabla_{\mathbf{r}}^{2}$ is the kinetic energy operator, $I$ is the index for atoms, $V_{I}^{local}$ and $V_{I}^{KB}$ are the local and Kleinman-Bylander parts of the pseduopotential for the $I$th atom, $V^{H}$ is the Hartree potential, $V^{XC}$ is the exchange-correlation (XC) potential and $V^{ext}$ is the potential of external field. Details in the calculation of $\braket{\xi_{i\alpha}(\mathbf{r}+\mathbf{R}_{s}+\mathbf{b}_{i})}{\hat{H}}{\xi_{j\beta}(\mathbf{r}+\mathbf{b}_{j})}$ are described in Ref. Ordejón, Artacho, and Soler, 1996. Within adiabatic local density approximation (LDA) and generalized gradient approximations (GGA) Yabana and Bertsch (1996) for the exchange-correlation functional, $V^{XC}$ does not depend explicitly on time $t$, i.e. $V^{XC}[\rho(\mathbf{r},t),t]=V^{XC}[\rho(\mathbf{r},t)]$. Thus, most XC functionals in ground-state DFT such as Perdew-Wang Perdew and Zunger (1981), Perdew-Burke-Ernzerhof Perdew, Burke, and Ernzerhof (1996), Becke-Lee-Yang-Parr Becke (1988); Lee, Yang, and Parr (1988), and van der Waals density functional Dion et al. (2004); Román-Pérez and Soler (2009) are compatible in this implementation of rt-TDDFT.

II.2 External field

To simulate the laser-matter interactions, time-dependent electric field $\mathbf{E}(t)$ is introduced to the Hamiltonian to represent the external time-dependent laser field in two different scenarios: the length gauge and vector gauge.

Within the length gauge, the effect of electric field $\mathbf{E}(t)$ is added to $V^{ext}$ as a scalar potential

[TABLE]

Time dependent $\mathbf{E}(t)$ can be tuned adopting any shape in time evolution. A most popular example is using the shape of a Gaussian wave packet

[TABLE]

where $f$ is the laser frequency, $t_{0}$ is the peak time, and $\phi$ is the phase factor.

We note that, the translational symmetry of Hamiltonian is broken by the introduction of finite external field $\mathbf{E}$ in the length gauge, since

[TABLE]

Thus, a common solution is using a sawtooth field along spatial direction ${\mu}\in\{x,y,z\}$

[TABLE]

where $L_{\mu}$ is the length of unit cell along ${\mu}$ and $\epsilon\rightarrow 0$. Thus, $-{E_{\mu}(t)}L_{\mu}/2\epsilon\rightarrow\infty$, which requires that charge density vanishes $\rho(x_{\mu})=0$ in the region $-\epsilon<x_{\mu}<+\epsilon$, otherwise the energy diverges. Thus, a vacuum layer is essential along ${\mu}$. The requirement for a vacuum layer limits the application of theoretical approaches using the length gauge field to study the extended systems. Since there is no vacuum layer in the extended bulk systems, the translational symmetry of the Hamiltonian is broken, $H(\mathbf{r}+\mathbf{R}_{s})\neq H(\mathbf{r})$, using the length gauge field. Plus, length gauge field is invalid in large systems and in short wavelength perturbation Lestrange, Egidi, and Li (2015).

Dynamical electric field in the vector gauge by introducting vector potential $\mathbf{A}$ could preserve the transitional symmetry of Hamiltonian, thus removes the requirement of the vacuum layer. Yabana et al. (2006, 2012) The relation between $\mathbf{E}$ and $\mathbf{A}$ is

[TABLE]

The Hamiltonian with the presence of $\mathbf{A}$ is then

[TABLE]

where

[TABLE]

within Rydberg atomic unit, where $e=\sqrt{2}$, $\hbar=1$ and $t=\hbar$/Ry. The unit of $\mathbf{k_{A}}$ is $\mathrm{Bohr}^{-1}$, the same as the unit of $\mathbf{k}$.

II.3 Propagation

With time-dependent (TD) Hamiltonian and overlap matrix, TDKS equation is solved to obtain $\ket{u_{n\mathbf{k}}(\mathbf{r},t+\Delta t)}$ from the state $\ket{u_{n\mathbf{k}}(\mathbf{r},t)}$ at the previous time step:

[TABLE]

where $\Delta t=t_{2}-t_{1}$ is the length of time step, $\ket{u_{n\mathbf{k}}(\mathbf{r},t)}$ is Bloch function and $t^{\prime}\approx({t_{1}+t_{2}})/{2}$.

It is rather difficult to evaluate $H_{\mathbf{k}}(t^{\prime})$ and $S_{\mathbf{k}}(t^{\prime})$ directly. Because $\Delta t$ is quite small ($<0.05$ fs), the ion positions barely changes from $t_{1}$ to $t_{2}$. Since $S_{\mathbf{k}}(t)$ is only determined by ionic positions (Eq. (1)), it is accurate enough to assume $S_{\mathbf{k}}(t^{\prime})\approx S_{\mathbf{k}}(t_{2})$. However, $H_{\mathbf{k}}(t)$ may largely change due to the rapid evolution of electrons. To approximate $H_{\mathbf{k}}(t^{\prime})$ properly, mid-point technique has been widely used Li et al. (2005a); Goings, Lestrange, and Li (2018).

Note that, $\ket{u_{n\mathbf{k}}(\mathbf{r})(t_{2})}$ is not explicitly dependent on other TDKS orbitals $\ket{u_{n^{\prime}\mathbf{k}^{\prime}}(\mathbf{r})(t_{1})}$ ($n^{\prime}\neq n$ or $\mathbf{k}^{\prime}\neq\mathbf{k}$), as a result of the $v$-representativity of the TDKS equations Runge and Gross (1984); Marques et al. (2012). It decouples the evolution equations of different TDKS orbitals and make TDDFT calculations practical. However, it nevertheless can account for both interband and intraband scatterings. Because $H_{\mathbf{k}}$ is determined by the total charge density, which is a weighted summation of all the occupied orbitals, there still exists an implicit coupling between different TDKS orbitals.

Numerically, the time propagator $\exp(-iS^{-1}_{\mathbf{k}}H_{\mathbf{k}}\Delta t)$ in Eq. (11) is expanded using first-order Crank-Nicholson scheme:

[TABLE]

Technically, since computing $S^{-1}_{\mathbf{k}}$ is the most time-consuming part in the calculation of Eq. (12), we minimize the times for its computing. $S^{-1}_{\mathbf{k}}$ is only updated when atomic positions, thus the center of NAOs, $\mathbf{b}_{i}$ are changed. Consequently, when ions are fixed, $S^{-1}_{\mathbf{k}}$ is computed only once at the first ionic step. Even with ions moving, $S^{-1}_{\mathbf{k}}$ only need to be updated once for each ionic step.

II.4 Updating charge density

With $\ket{u_{n\mathbf{k}}(\mathbf{r},t_{2})}$ solved in Eq. (11), the density matrix (DM) $\rho_{i\alpha,j\beta}(t_{2})$ is computed accordingly as:

[TABLE]

where $q_{n,\mathbf{k}}$ is the electronic population of the band $n$ at $\mathbf{k}$, and $c_{n,j\beta,\mathbf{k}}(t_{2})$ is the coefficient of $\ket{u_{n\mathbf{k}}(\mathbf{r},t_{2})}$ in NAO basis:

[TABLE]

II.5 Self-consistent evolution

We use the self-consistent process described in Ref. Meng and Kaxiras (2008) during the time evolution of charge density. This process substantially increases the numerical stability Ren, Kaxiras, and Meng (2010). All the criteria for convergence test developed in SIESTA are compatible with the current approach, such as using the maximum element of the DM difference, the energy difference, or the Harris energy difference, etc. as a criterion for achieving self-consistency. Soler et al. (2002).

Here, we use DM difference as an example. Convergence in charge density during time evolution is reached when

[TABLE]

where $\eta$ is about 10*-4*.

II.6 Mixing

If not converged, the linear mixing of DM is needed to generate the input DM for computing charge density $\rho_{next}$ at the next cycle, instead of using $\rho_{new}$ directly,

[TABLE]

where the $\rho$ on the right side of Eq. (16) is the input DM and $\rho_{new}$ is the output DM, and $w$ is the mixing weight, usually $w=0.1-0.5$.

II.7 Postprocessing

If self-consistent time evolution of charge density is converged, the postprocessing steps are evoked, including the calculation of total energy, Hellmann-Feynman forces, ionic movements, etc. These functions are implemented in SIESTA Soler et al. (2002) and compatibly used in TDAP Meng and Kaxiras (2008). We note that, rt-TDDFT in atomic orbital basis gives rise to additional Pulay terms that contribute to the force evaluationsDing et al. (2015b); Isborn2007a; Schlegel et al. (2001). The total force is the combination of Hellmann-Feynman force and Pulay term. With the calculated forces, the coupled electron-ion motion can be simulated based on classical ionic trajectories, in the framework of Ehrenfest dynamics. In Ehrenfest dynamics, the forces on the ions are averaged over the adiabatic electronic states along all possible ionic paths. If one path is dominating or many similar potential energy surfaces are involved, Ehrenfest dynamics works very well Tully (1990); otherwise, classical trajectory approximations in Ehrenfest dynamics become less accurate Parandekar and Tully (2006); Hack and Truhlar (2000). Furthermore, detailed balance for quantum electronic states is absent in the Ehrenfest dynamics. Thus, the applications of the present method are limited to the cases where the averaged potential energy surfaces yields a reasonable description of coupled electron-ion dynamics. Since we focus on the dynamics of excited electrons in this work, the ions are fixed in the simulations.

Here we introduce some analysis in detail for typical rt-TDDFT simulations. First, we could evaluate the state-to-state transition probabilities between TDKS orbitals during time evolution Li et al. (2005a); Rohringer, Peter, and Burgdörfer (2006):

[TABLE]

where $\ket{v_{n\mathbf{k}}}$ is the adiabatic basis satisfying

[TABLE]

The population $\mathcal{q}_{n\mathbf{k}}$ of the adiabatic state $n\mathbf{k}$ is thus projected from the TDKS orbitals at a given time as:

[TABLE]

where $n_{\mathbf{k},occ}$ is the occupied state at $\mathbf{k}$ point.

For finite systems and surface slabs, we can calculate time-dependent dipole moment along the direction. For periodic systems, the dipole moment is ill-defined. Instead, we calculate time dependent current,

[TABLE]

as the response function.

III Results and Discussion

III.1 Momentum-resolved versus supercell approaches

To demonstrate the $\mathbf{k}$-resolved algorithm, we choose graphene as the model system (see Fig. 2(a)). An exotic property of graphene (also of other Dirac materials) is the linear dispersion near K point, namely, $E(\mathbf{k})=v_{F}\mathbf{k}$, where $E$ the band energy and $v_{F}$ is the Fermi velocity which could reach $10^{6}$ m/s. To describe all the Bloch electrons, especially those near the Fermi energy, two kinds of strategies are used: unit cell calculations with $\mathbf{k}$-resolved reciprocal space samplings or a supercell approach with $\Gamma$-only $\mathbf{k}$-sampling. To demonstrate the advantages of the $\mathbf{k}$-resolved algorithm, we compare three cases:

(i) unit cell with the Monkhorst-Pack Monkhorst and Pack (1976) $N_{k}\times N_{k}\times 1$ $\mathbf{k}$ point mesh, to cover all important special $\mathbf{k}$-points $M$, $\Gamma$ and $K$, facilitating a line-mode analysis along $M\rightarrow\Gamma\rightarrow K\rightarrow M$ [Fig. 2(b)];

(ii) $N_{c}\times N_{c}\times 1$ supercell with single $\Gamma$ point; and

(iii) $N_{c}\times N_{c}\times 1$ supercell with single $K$ point.

To compare the computation accuracy of these three cases, we define an error function $\Delta$ as,

[TABLE]

where $T$ is the total simulation time, $E^{ref}_{ex}$ is the excitation energy of the reference case and $E_{ex}(t)$ is the excitation energy

[TABLE]

where $E_{KS}$ is the total energy of the system.

Here, we evaluate $\Delta$ under such settings: the Gaussian-shaped laser pulse [Eq. (5)] with $\phi=0$, $t_{0}=7.0$ fs, $\sigma=2.0$ fs, $f=21.93$ eV is applied; the total simulation time is $T=20$ fs; and the reference energy $E^{ref}_{ex}$ is calculated with $60\times 60\times 1$ $\mathbf{k}$-point mesh. A diagram to illustrate the definition of $\Delta$ is shown in the inset of Fig. 3. The time step is chosen to $\Delta t=0.02$ fs and the total time is $20$ fs. Troullier-Martin pseudopotentials Troullier and Martins (1991), adiabatic local density approximation (ALDA) exchange-correlation functional Perdew and Zunger (1981); Yabana and Bertsch (1996) and an auxiliary real-space grid equivalent to a plane-wave cutoff of $75$ Ry are used. In description of C atoms, we use a basis set of 8 double-$\zeta$ orbitals {2s(2$\zeta$), 2px(2$\zeta$), 2py(2$\zeta$), 2pz(2$\zeta$)} and 5 polarization orbitals { $\mathrm{P_{d_{xy}}}$, $\mathrm{P_{d_{yz}}}$, $\mathrm{P_{d_{z^{2}}}}$, $\mathrm{P_{d_{xz}}}$, $\mathrm{P_{d_{x^{2}-y^{2}}}}$ }. We calculate the test cases with one 8-core Intel(R) Xeon(R) CPU E5-2650@2.00GHz.

We plot $\Delta$ of these three cases in Fig. 3. The error $\Delta$ decreases as $N=N_{k}$ (or $N_{c}$) increases. The absolute value of $\Delta$ on the the same scale is achieved with $N_{c}=N_{k}$. That is to say, the unit cell approach with $N\times N\times 1$ $\mathbf{k}$-point mesh is as accurate as the approach using a $N\times N\times 1$ supercell. To achieve an accuracy with the $\Delta\leq 2$ meV/atom, $N_{k}=24$ is needed. Thus, it can be predicted that $N_{C}=24$ is needed for the supercell approach.

However, we emphasize that the computational cost for calculating $N_{c}\times N_{c}\times 1$ supercell is extremely heavy. As shown in Fig. 4, solving Eq. (11) dominates ($\sim 80$%) the computer time consumption at large $N_{k}$ ($N_{c}$), which scales linearly with the total number of $\mathbf{k}$ points, $N_{k}^{2}$, and quadratically with the total number of atoms, $N_{c}^{2}$. The CPU clock time $t_{c}$ approximately scales as O($N^{2}_{k}\times N_{c}^{4}$). Thus, $t_{c}=N_{c}^{4}$ for supercell calculation, while $t_{c}=N_{k}^{2}$ for$\mathbf{k}$-resolved calculations at the same level of accuracy.

As $N$ increases, this difference become more significant. For supercell calculations, we are able to only compute supercells up to $N_{C}=9$, which already costs over $2\times 10^{3}$ min. At the same accuracy level, $N_{k}=9$ calculation costs only $20$ min, which is only 1/100 of that for $N_{C}=9$ case, consistent with the time complexity analysis $N_{k}^{2}/N_{c}^{4}=1/81$. As mentioned above, $N_{k}=24$ or $N_{C}=24$ is needed for relatively accurate calculations. To fulfill this requirement, calculation with $N_{k}=24$ costs only about 1 hour, showing that it is readily accessible and efficient. In contrast, calculating a $N_{c}=24$ supercell would require a computer time over 576 hours (24 days) and thus heavy in real applications. Regarding the computational accuracy and efficiency, we choose $\mathbf{k}$-point mesh $24\times 24\times 1$ to achieve an extremely dense sampling of the Brillouin zone.

With the small unit cell of graphene, the number of real space grids is $\sim$1000, which is 30 times of the number of NAOs used. Considering the evolution algorithm has the computational complexity of $O(n^{2})$, where $n$ is the number of basis functions, the computer time for wavefunction evolution using NAO basis is largely reduced to 1/90 of that using real space grid basis. In practical calculations using the same number of message-passing-interface (MPI) processes, the reduction in the total computer time is tested to be about 1/5 to 1/10 depending on the systems under consideration Meng and Kaxiras (2008); Lian et al. (2018).

III.2 Out-of-plane excitation

We then adopt a laser field perpendicularly polarized to the graphene plane to excite electrons in graphene, i.e. in a set-up of small angle scattering. Since there is a vaccum layer along the out-of-plane direction, the laser field in the length gauge can be used.

We first calculate the dielectric function of graphene to locate the photon energy for resonant excitation, $\alpha_{\mu,\nu}$, where $\mu$, $\nu$ denote the spatial direction $\mu$, $\nu$ $\in$ {$x$, $y$, $z$ }. The $\alpha_{\mu,\nu}$ describes the response of dipole moment $P_{\mu}(\omega)$ to the electric field $E_{\nu}(\omega)$ in the frequency domain,

[TABLE]

In rt-TDDFT calculations, we apply the electric field $E_{\nu}(t)$ and obtain the dipole moment $P_{\mu}(t)$ in time domain. Then we carry out the Fourier transform to get Eq. [23],

[TABLE]

We then obtain

[TABLE]

In principle $E_{\nu}(t)$ can be in an arbitrary shape with time. However, in practice, it is better to choose Dirac function $E^{\delta}_{\nu}(t)=E_{\nu 0}\delta(t)$, or the Heaviside step function $E^{\theta}_{\nu}(t)=E_{\nu 0}[1-\theta(t)]$ to include components $E_{\nu}(\omega)$ at all $\omega$, since we have

[TABLE]

Here we choose the latter form:

[TABLE]

which leads to

[TABLE]

Importantly, $\mathrm{Im}\{\alpha_{\mu,\mu}(\omega)\}$ characterizes the optical absorbance at $\omega$ along the $\mu$ direction.

We calculate the imaginary part of the dielectric function along the out-of-plane $z$ direction of graphene, $\mathrm{Im}\{\alpha_{z,z}(\omega)\}$. As shown in Fig. 5, $\mathrm{Im}\{\alpha_{z,z}(\omega)\}$ are almost the same with the increase of $E^{0}$ from $0.05$ to $0.5$ V/Å, indicating the linear response theory is appropriate in this range of light illumination. The absorption peaks are located at relative high energies ($>20$ eV). The first absorption peak is located at $21.93$ eV. We choose this photon energy to simulate the resonant excitation of graphene in the perpendicular direction.

We demonstrate the excitation dynamics of graphene at a resonant light frequency of $\omega_{r}=21.93$ eV, and compare to the case at the non-resonant light frequency of $\omega_{nr}=2.0$ eV. We characterize the overall excitation through tracking the number of excited electrons, as well as the total energy change during the excitation process as a function of time. The number of excited electrons $n(t)$ is calculated as,

[TABLE]

where $\mathcal{q}_{n\mathbf{k}}(t)$ is obtained from Eq. (19), and ${unocc}$ denotes the unoccupied TDKS states.

As shown in Fig. 6, different behaviors are observed for the two excitation conditions. The excited electrons $n(t)$ and excitation energy $E_{ex}(t)$ increases at $\omega_{r}$, while no response is observed at $\omega_{nr}$. The same results are obtained at other non-resonant light frequencies of $1.0$, $2.0$, $4.0$ eV. It verifies that the calculated $\mathrm{Im}\{\alpha_{z,z}(\omega)\}$ characterizes well selectivity in optical absorption: only the light with the right photon energy $\omega$, at which $\mathrm{Im}\{\alpha_{z,z}(\omega)\}$ peaks, has a strong absorption.

We discuss the resonant case here. In general, $n(t)$ is similar to the shape of the laser pulse, while two special features are observed. Firstly, the time variation in $n(t)$ has a $1.4$ fs delay from the laser field. This delay represents the intrinsic response time of graphene to laser field, namely, the time needed for light absorption and electronic transitions. Secondly, $n(t)$ decreases but does not vanish after the end of light pulse. Thus, we propose that two kinds of excitation process exist: one is the transient excited electrons, which quickly vanishes after the laser pulse is off; another is the residual excited electrons, which live relatively longer. Residual $n(t)$ would decrease with the occurrence of electron-electron and further electron-phonon scatterings at the time scale of $100$ fs, thus is not observed in our short-time simulation ($<20$ fs). We note that, the dependence on history is absent in the calculations with adiabatic exchange-correlation functionals, which causes less accurate prediction of the lifetime of excited states and ionic forces on a long time scale Maitra, Burke, and Woodward (2002); Elliott et al. (2012); Maitra ; Ullrich (2006); Agostini et al. (2015).

To verify our assumption about the existence of two kinds of excitation processes, we further distinguish the excitation with $\mathbf{k}$-point resolution. We choose six snapshots of $\mathcal{q}_{n\mathbf{k}}(t)$ defined in Eq. (19), as shown in Fig. 7. At $t=2.0$ fs with the absence of laser pulse, no excitation is observed at all $\mathbf{k}$ points. At $t=4.0$ fs, the excitation is still ignorable, while the laser field is just turned on, due to the delay in electronic response we discussed above. At the peak time of the laser pulse $t=7.0$ fs, $\mathcal{q}_{n\mathbf{k}}(t)$ shows a significant distribution over many $\mathbf{k}$-points. We mark the dominant excitation mode as L, which involving bonding $\pi$ and antibonding $\pi$ bands. With $t$ increases from $7$ fs to $12$ fs, the L mode excitation rapidly decreases. In contrast, two new modes (labeled by their locations in the reciprocal space, K1 and K2) increases and become dominant. K1 and K2 modes maintain within $20$ fs while L mode gradually vanishes. Thus, with the assistance of newly developed $\mathbf{k}$-resolved algorithm, we are able to successfully distinguish these two kinds of excitation processes: L mode produce the transient excited electrons while K1 and K2 modes produce the residual excited electrons.

Although K1 and K2 are both long-lived excitations, their time-dependence is quite different. We plot $\mathcal{q}_{n\mathbf{k}}(t)$ as a function of $t$ at three $\mathbf{k}$ points $\Gamma$, K1 and K2, as shown in Fig. 8. For L mode (represented by photoexcitation at $\Gamma$ point), the clear transient character is demonstrated. The excitation only exists when the laser field is present, consistent with the observations in Fig. 7. However, for K1 and K2 modes, new differences are observed. Excited electrons in K1 mode increases monotonically, while $\mathcal{q}_{n\mathbf{k}}(t)$ at K2 shows an oscillation with a periodicity of $T_{K_{2}}\sim 5$ fs. These different behaviors are due to different excitation energies of three modes, originated from different band structures at the different $\mathbf{k}$-point. For instance, the oscillation of $K_{2}$ mode is analogous to the beating,

[TABLE]

For K2 mode, $\omega_{K_{2}}=22.75$ eV is the energy difference between the two electronic bands involved in the optical transition at K2 (initial and final states), and $\omega_{2}=\omega_{r}=21.93$ eV is the driving photon energy. Beat frequency $T_{b}={4\pi}/{(\omega_{1}-\omega_{2})}=5.07$ fs, which is close to the observed oscillation periodicity $T_{K_{2}}$. Thus, the oscillation of K2 mode is the beat formed by the intrinsic band energy difference and the driving laser frequency. In contrast, $K_{1}$ mode excitation has very close energies: $\omega_{K_{1}}=21.59$ eV and $\omega_{r}=21.93$ eV, thus only a half period of the beat ($T_{K_{1}}=12.4$ fs) is observed in our simulation. For L mode, the excitation energy is $19.32$ eV, far below the $\omega_{r}$. A non-resonant interference shows up instead of beating. The rich photoexcitation phenomena discussed above and the associated complex dynamic behaviors hint for the needs for developing efficient rt-TDDFT algorithms with momentum resolution. By introducing a new degree of freedom in the reciprocal space, the $\mathbf{k}$-resolved dynamics labels the distinct excitation processes as well as final distribution of excited states in the Brillouin zone after the incidence of laser pulses.

III.3 In-plane excitation

For a laser pulse with its field polarization lying parallel to the atomic plane of graphene (i.e., normal incidence), adding a vacuum layer along the laser polarization direction is not possible for the periodical extended system such as graphene. Therefore we adopt the vector potential approach to simulate the in-plane laser-graphene interaction. The graphene sheet is illuminated with a linearly polarized laser pulse, as shown in Fig. 9. We note that in-plane excitation is well described by the Fermi’s golden rule. Only the bands with an energy gap $\Delta E_{g}(\mathbf{k})$ equal to the photon energy $\omega$ will be excited. As a result, in-plane polarized laser excites electrons near the Dirac point for photon energies $\leq$5 eV, see Fig. 9. The momentum-resolved simulation will distinguish the photoexcitation induced by a laser pulse with different photon energies $\omega$.

Here, we use four different wavelengths $\lambda=$ 1200 nm, 600 nm, 400 nm, and 300 nm for the laser pulse, corresponding to photon energies $\omega=$ 1.03, 2.06, 3.10, 4.13 eV, respectively, to excite graphene in the in-plane direction. For simplicity the laser field is polarized perpendicular to the C-C bond of graphene lattice (referred to as $y$ direction). The momentum resolved excitation patterns in the reciprocal space are shown in Fig. 10 (a), with the corresponding band energy difference $\Delta E_{g}(\mathbf{k})$ shown in Fig. 10 (b). It is clear that only the $\mathbf{k}$ points with $\Delta E_{g}(\mathbf{k})=\omega$ are excited. This agreement justifies the validation of the vector gauge used in the current TDDFT implementation.

Furthermore, we note that the presence of strong laser field breaks the six-fold rotational symmetry of the graphene lattice. For instance, with $\omega=4.13$ eV, photoexcitation at two $M^{\prime}$ points are absent, while excitations at other symmetric $M$ points are observed. This symmetry breaking is caused by presence of linearly polarized laser field along the $y$ direction.

It can be explained with a two band model of graphene (see appendix). The excellent agreement on the excitation outcome between the model Hamiltonian and first-principles quantum dynamics simulations justify the validity of our rt-TDDFT algorithm with a vector gauge field. We therefore expect that it is readily applicable to investigate quantum dynamics of a variety of electronic phases such as charge/spin density waves, Mott insulators, valley electronics, and electronic melting in two-dimensional materials and conventional semiconductors.

To demonstrate the general applicability of the present approach, we tackle photoexcitation induced electron dynamics in a complex material. The layered transition-metal dichalcogenides such as 1T-TaS2 have been widely studied in literature to understand charge density wave (CDW) physics in real materials, whose structure is shown in Fig. 11(a). The 1T-TaS2 is a typical quasi two-dimensional CDW material with a pristine lattice constant of 3.36 Å in the undistorted 1T phase. In ground state, the lattice undergoes a structural reconstruction forming a $\sqrt{13}\times\sqrt{13}$ superstructure with star-of-David patterns. Laser induced phase dynamics in 1T-TaS2 has been investigated in recent experiments, where its responses to ultrashort laser pulses play a critical role. Here we study the carrier distribution in 1T-TaS2 upon ultrafast laser excitation. As shown in Fig. 11(b), the excitation energy strongly oscillates with the field of laser pulse. The excitation energy deposited by the laser pulse is $\sim$12 eV/cell after laser illumination with a photon energy of $\hbar\omega=1.55$ eV and pulse width of 8 fs. The carrier distribution at 20 fs after the passing of the laser pulse is shown in Fig. 11(c). The majority of excited electrons and holes are located at energies ranging from $-$2 to 2 eV near the Fermi level. It indicates that the photoexcitation mainly consists of single-photon processes as well as a minor fraction of two-photon processes (with excited electrons located at $\sim$3 eV and holes at $-3$ eV).

IV Conclusions

In conclusion, we have developed $\mathbf{k}$-resolved rt-TDDFT algorithms using efficient numerical atomic basis. It enables large-scale rt-TDDFT simulations of extended systems including solids, interfaces, and two-dimensional materials with a rather small unit cell, significantly reducing the heavy computational cost of typically rt-TDDFT simulations. Consequently, $\mathbf{k}$-resolved excitation dynamics in periodical crystal materials are observed. The key advantages of this unique approach includes:

i) The $\mathbf{k}$-resolved real-time evolution algorithm introduces the important $\mathbf{k}$-space resolution and a new degree of freedom, which is essential to describe key quantities and important physics in photoexcited condensed matter materials. The use of many $\mathbf{k}$-points with a rather small unit cell also significantly improves the computational efficiency of rt-TDDFT calculations for photoexcitation in solids.

ii) Different from approaches using real space grids and all-electron full-potential linearized augmented-planewaves, the adoption of numerical atomic basis in the present implementation reduces the number of required basis functions to one-hundredth of its original value, making rt-TDDFT computation of realistic large systems (comprising $\sim$500 atoms and lasting for $\sim$1000 fs) plausible. In addition, with a relatively small real-space cutoff for NAOs, the order-$N$ linear scaling with respect to the system size can be achieved.

iii) Both electronic and ionic degree of freedoms are evolved, therefore a complete information on electronic wavefunctions and ionic movements during real time evolution can be provided for simulations of complex materials and rich phenomena far from equilibrium.

When applied to study photoexcitation dynamics of a prototypical model material–graphene, the $\mathbf{k}$-resolved algorithm enables the observation of $\mathbf{k}$ selective excitation modes. Three distinct modes are excited, located at different $\mathbf{k}$. In-plance excitation of the Dirac electrons in graphene can be understood by assuming an effective vector field of laser field, via taking into account the angular dependence of optical transition matrix elements. This kind of $\mathbf{k}$ dependent electronic dynamics are ubiquitous in solids. Thus, $\mathbf{k}$-resolved rt-TDDFT algorithm is an important development for investigating ultrafast photoexcitation dynamics and electron-electron scattering, and is expected to be widely used in the future.

V acknowledgement

We acknowledge partial financial support from MOST (Grant Nos. 2016YFA0300902 and 2015CB921001), NSFC (Grant Nos. 11774396, 11474328, and 91850120), and CAS (Grant No. XDB07030100).

VI APPENDIX: THE TWO-BAND MODEL OF GRAPHENE

The ground state Hamiltonian of two-band model of graphene reads,

[TABLE]

where $k_{x}$, $k_{y}$ is the $\mathbf{k}$ coordinate, $\sigma$ is the Pauli matrices, $v_{F}$ is the Fermi velocity. $v_{F}=1$ eV$\cdot$Bohr for simplification. The units of $k_{x}$ and $k_{y}$ are chosen as Bohr*-1*. The energy unit is thus eV. The eigenvalues and eigenvectors are solved as,

[TABLE]

[TABLE]

Thus, the initial state wavefunction is the ground state $\psi(t=0)=\phi_{0}$.

A vector field polarized along $y$ can be introduced as,

[TABLE]

where $A(t)$ is the vector gauge field. The time-dependent Hamiltonian is thus,

[TABLE]

The wavefunction at time $t$ can be obtained from time-dependent Schrödinger equation,

[TABLE]

which can be expanded with $\ket{\phi_{0}}$ and $\ket{\phi_{1}}$ basis,

[TABLE]

where

[TABLE]

is the time-dependent coefficients. All equations are solved numerically with the qutip package Johansson, Nation, and Nori (2012, 2013).

We can reproduce the symmetry breaking in the distribution of excited electrons in $\mathbf{k}$ space induced by linearly polarized laser. We analyze the coefficients of $|c_{2}(t)|^{2}$ with $k_{x}=\cos\theta,k_{y}=\sin\theta$, under the vector field $\mathbf{A}=0.2$ Bohr*-1*, where $\theta$ is the angle between $\mathbf{k}$ and $\mathbf{A}$, as shown in Fig. 12. Thus, the energy difference $\Delta E_{g}(\mathbf{k})=2.0$ eV. These $\mathbf{k}$ points are only excited with $\omega=2.0$ eV, consistent with the results from TDDFT and Fermi’s golden rule.

To explain the origin of the symmetry breaking, the excited electrons at different $\mathbf{k}$ points at the end of laser pulse $|c_{2}(t=50~{}\mathrm{fs})|^{2}$ are shown in Fig. 12(a). It suggests that, the effect of linearly polarized laser on point $\mathbf{k}$ is not solely characterized by the $\mathbf{A}$ field, but also related to the angle $\theta$ between $\mathbf{k}$ and $\mathbf{A}$. With $\theta=0$ and $\pi$, i.e. the $\mathbf{k}$ is parallel/anti-parallel to the $\mathbf{A}$ field, the excitation is fully suppressed, while the excitation is the maximum with $\theta=\pi/2$ and $3\pi/2$. An effective field $A_{\mathrm{eff}}=A\sin\theta$, always perpendicular to the vector $\mathbf{k}$, is thus introduced to induce electronic transitions at $\mathbf{k}=(k_{x}=\cos\theta,k_{y}=\sin\theta)$, as shown in Fig. 12(b). It explains the origin of the symmetry breaking in TDDFT simultions (Fig. 10). The excitations at $\mathbf{k}$ points are the results of the combined effects of energy match and the angle $\theta$ between the $\mathbf{k}-\mathbf{K}$ and $\mathbf{A}$ field, where $\mathbf{K}$ is the coordinates of the adjacent Dirac point. Since $\mathbf{A}$ field is along $y$, $\sin\theta=0$ for all the $\mathbf{k}$ points with $\mathbf{k}-\mathbf{K}$ parallel to the polarization direction. Thus, this is no effective field to introduce photoexcitations at the two $M^{\prime}$ points.