Ultrafast nonequilibrium evolution of excitonic modes in semiconductors

TL;DR

This paper investigates the ultrafast nonequilibrium dynamics of excitonic modes in semiconductors using advanced numerical methods, revealing distinct behaviors based on excitation types and demonstrating excitonic condensation out of equilibrium.

Contribution

It introduces a comprehensive comparison of numerical methods for simulating excitonic dynamics and uncovers new phenomena related to collective modes and nonequilibrium excitonic condensation.

Findings

GKBA aligns best with iTEBD in simulations.

Different excitations lead to distinct excitonic behaviors.

Evidence of excitonic condensation out of equilibrium.

Abstract

We study the time evolution of excitonic states after photo-excitation in the one-dimensional spin-less extended Falicov-Kimball model. Several numerical methods are employed and benchmarked against each other: time-dependent mean-field simulations, the second-Born approximation (2BA) within the Kadanoff-Baym formalism, the generalized Kadanoff-Baym Ansatz (GKBA) implemented with the 2BA and the infinite time-evolving block decimation (iTEBD) method. It is found that the GKBA gives the best agreement with iTEBD and captures the relevant physics. We find that excitations to the particle-hole continuum and resonant excitations of the equilibrium exciton result in a qualitatively different dynamics. In the former case, the exciton binding energy remains positive and the frequency of the corresponding coherent oscillations is smaller than the band gap. On the other hand, resonant excitations trigger a collective mode whose frequency is larger than the band gap. We discuss the origin of these different behaviors by evaluating the nonequilibrium susceptibility using the nonthermal distribution and a random phase approximation. The peculiar mode with frequency larger than the band gap is associated with a partial population inversion with a sharp energy cutoff. We also discuss the effects of the cooling by a phonon bath. We demonstrate the real-time development of coherence in the polarization, which indicates excitonic condensation out of equilibrium.

Equations

$$\hat{H}(t)$$

$${\hat{H}}_{\mathrm{kin}}$$

$${\hat{H}}_{\mathrm{int}}$$

$$\displaystyle J_{a}({\bf r}_{ij},t)=J_{a}({\bf r}_{ij})\exp\Big{[}i\frac{q}{\hbar}{\bf r}_{ij}\cdot{\bf A}(t)\Big{]},$$

$${\hat{H}}_{\mathrm{dip}}(t)=-\mathbf{E}(t)\cdot \sum _i{\hat{\mathbf{P}}}_i=-\mathbf{E}(t)\cdot \sum _{i,a}{\mathbf{d}}_a{\hat{c}}_{i,a}^{†}{\hat{c}}_{i,\bar{a}}.$$

$${\hat{H}}_{\mathrm{kin}}(t)$$

$${\hat{H}}_{\mathrm{dip}}(t)$$

$${\mathbf{h}}_{\mathrm{kin},\mathbf{k}}(t)$$

$${\mathbf{h}}_{\mathrm{dip},\mathbf{k}}(t)$$

$${\rho }_{ia,jb}(t)$$

$${\rho }_{\mathbf{k},a,b}(t)$$

$${\chi }_{BA}^R(t,t^{\prime })=\sum _{ijlm}{B}_{ml}{\chi }_{ml,ij}^R(t,t^{\prime })A_{ij},$$

$${\chi }_{ml,ij}(t,t^{\prime })=$$

$$+i⟨{\mathcal{T}}_{\mathcal{C}}{\hat{c}}_l^{†}(t){\hat{c}}_m(t)⟩⟨{\mathcal{T}}_{\mathcal{C}}{\hat{c}}_j^{†}(t^{\prime }){\hat{c}}_i(t^{\prime })⟩$$

$${\hat{H}}_{\mathrm{MF}}(t)=\sum _{\mathbf{k}}\left[\begin{array}{cc}{\hat{c}}_{\mathbf{k},c}^{†}& {\hat{c}}_{\mathbf{k},v}^{†}\end{array}\right]\cdot {\mathbf{h}}_{\mathrm{MF},\mathbf{k}}(t)\cdot \left[\begin{array}{c}{\hat{c}}_{\mathbf{k},c}\\ {\hat{c}}_{\mathbf{k},v}\end{array}\right],$$

$${\mathbf{h}}_{\mathrm{MF},\mathbf{k}}(t)$$

$${\mathbf{h}}_{\mathrm{Hartree},\mathbf{k}}(t)$$

$${\mathbf{h}}_{\mathrm{Fock},\mathbf{k}}$$

$${\mathbf{\chi }}^R(\omega ;\mathbf{q})=[1-{\mathbf{\chi }}_0^R(\omega ;\mathbf{q})\mathbf{\Theta }{]}^{-1}{\mathbf{\chi }}_0^R(\omega ;\mathbf{q}).$$

$${\chi }_{0,\mu \nu }(t;\mathbf{q})=-i\theta (t)\frac{1}{N}$$

$$\displaystyle-{\rm tr}[\boldsymbol{\sigma}_{\mu}\boldsymbol{\mathcal{G}}^{<}_{\bf k+q}(t)\boldsymbol{\sigma}_{\nu}\boldsymbol{\mathcal{G}}^{>}_{\bf k}(-t)]\bigl{\}}.$$

$${\mathcal{G}}_{aa,\mathbf{k}}^{<}(t)$$

$${\mathcal{G}}_{aa,\mathbf{k}}^{>}(t)$$

$${\chi }_{0,\mu \nu }(\omega ;\mathbf{q})=\frac{1}{N}\sum _{\mathbf{k},a,b}\frac{\mathrm{tr}[{\mathbf{W}}_a{\mathbf{\sigma }}_{\mu }{\mathbf{W}}_b{\mathbf{\sigma }}_{\nu }](n_a(\mathbf{k}-\mathbf{q})-n_b(\mathbf{k}))}{\omega +i{0}^{+}-(E_b(\mathbf{k})-E_a(\mathbf{k}-\mathbf{q}))},$$

$$\mathbf{\gamma }={\mathbf{L}}^{-1}\mathbf{\chi }\mathbf{L},{\mathbf{\gamma }}_0={\mathbf{L}}^{-1}{\mathbf{\chi }}_0\mathbf{L},\mathbf{L}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& i\\ i& 1\end{array}\right].$$

$${\gamma }_{0,11}(\omega )$$

$${\gamma }_{0,22}(\omega )$$

$${\gamma }_{\mu \mu }(\omega )$$

$$G_{ab,\mathbf{k}}(t,t^{\prime })=-i⟨{\mathcal{T}}_{\mathcal{C}}{\hat{c}}_{\mathbf{k},a}(t){\hat{c}}_{\mathbf{k},b}^{†}(t^{\prime })⟩.$$

$$\left[i{\partial }_t-{\mathbf{h}}_{\mathrm{MF},\mathbf{k}}(t)\right]{\mathbf{G}}_{\mathbf{k}}(t,t^{\prime })$$

$${\mathbf{\Sigma }}^{2B}[\mathbf{G}](t,t^{\prime })={\mathbf{\Sigma }}^{2Bd}[\mathbf{G}](t,t^{\prime })+{\mathbf{\Sigma }}^{2Bx}[\mathbf{G}](t,t^{\prime }).$$

$${\Sigma }_{ab,\mathbf{k}}^{2Bd}(t,t^{\prime })$$

$$\times G_{\bar{b}\bar{a},\mathbf{p}}(t^{\prime },t),$$

$${\Sigma }_{ab,\mathbf{k}}^{2Bx}(t,t^{\prime })$$

$$\times G_{\bar{b}\bar{a},\mathbf{p}}(t^{\prime },t).$$

$${\chi }_{ml,ij}(t,t^{\prime })=-i\mathrm{tr}[{\mathbf{e}}_{lm}{\mathbf{G}}_0(t,t^{\prime }){\mathbf{e}}_{ji}{\mathbf{G}}_0(t^{\prime },t)]$$

Full text

I Introduction

Excitonic states play a central role in photo-excited semiconductors, nanostructures and molecules and have been studied extensively in the context of photo-voltaic applications Haug and Koch (1990); Ostroverkhova (2016); Scholes and Rumbles (2006); Koch et al. (2006) and charge migration.Hill et al. (2000); Falke et al. (2014); Boström et al. (2018) In particular, two-dimensional (2D) materials – especially transition metal chalcogenides (TMCs) – are currently attracting a lot of interest, fueled by the possibility of creating tailored heterostructures.Novoselov et al. (2005); Radisavljevic et al. (2011); Heine (2015); Novoselov et al. (2016) Due to the low dimensionality of TMCs, the Coulomb interaction is weakly screened, thus giving rise to pronounced interaction effects and excitonic features. TMCs exhibit large exciton binding energies, which can be of the order of a few hundred meV.He et al. (2014); Heine (2015); Cudazzo et al. (2016) Apart from the importance of excitons as excited states dominating the in-gap optical absorption – known as virtual or coherent excitons Haug and Koch (1990); Schäfer and Wegener (2013); Koch et al. (2006) – excitons can also be present in the ground state. For sufficiently large binding energy, these excitons can condense collectively, forming an excitonic insulator (EI).Kohn (1967); Jérome et al. (1967); Halperin and Rice (1968) Because of the strong Coulomb interaction, TMCs are among the best candidates for realizing the EI phase.Cercellier et al. (2007); Hellmann et al. (2012); Wakisaka et al. (2009); Kaneko et al. (2013); Lu et al. (2017)

While virtual excitons in semiconductors are usually considered in the linear response regime, stronger excitations and out-of-equilibrium dynamics have also been in the spot light. Dynamics of semiconductors after strong excitations and the realization of the EI phase there have been investigated theoretically,Comte and Mahler (1986); Östreich and Schönhammer (1993); Perfetto et al. (2019) and the relevant photo-dressed states have been observed recently.Murotani et al. (2019) Furthermore, the strong light-matter coupling in TMCs,Britnell et al. (2013) which can be enhanced by orders of magnitude in a micro-cavity setup,Peter et al. (2005); Bisht et al. (2019); Latini et al. (2019) implies that excitonic properties need to be investigated beyond linear response. Important examples for nonequilibrium setups also include the optical Stark effect Sie et al. (2015) and the ultrafast charge transfer in photo-excited bilayer TMCs.Ceballos et al. (2014); Hong et al. (2014) In addition, in order to unravel the mechanisms of the photo-induced enhancement Mor et al. (2017); Murakami et al. (2017) or ultrafast melting of EI orders Okazaki et al. (2018), it is essential to develop an understanding of the dynamics of bound electron-hole pairs in strongly photo-excited systems.Hellmann et al. (2012); Golež et al. (2016); Mor et al. (2017); Murakami et al. (2017); Tanabe et al. (2018); Tanaka et al. (2018); Okazaki et al. (2018)

In the linear response regime, excitons are typically treated within the framework of the Bethe-Salpeter equation (BSE) Rohlfing and Louie (1998); Albrecht et al. (1998); Cudazzo et al. (2016) in combination with the kernel determined by the Hartree-Fock self-energy (the random-phase approximation, RPA) or the $GW$ approximation. Extending the BSE to a nonequilibrium scenario is possible,Perfetto et al. (2015) but currently out of reach for realistic systems. Time-dependent approaches are a promising alternative route for computing the linear Berghäuser and Malic (2014); Perfetto et al. (2016); Murakami et al. (2016a, b); Sangalli et al. (2018) and beyond-linear response.Attaccalite et al. (2011, 2018) In particular, the nonequilibrium Green’s functions (NEGF) Stefanucci and Leeuwen (2013) approach provides a natural way of extending the many-body perturbation theory to the time domain. However, a priori it is unclear which scheme works best out of equilibrium. For instance, the spurious effects of fully self-consistent $GW$ Grumet et al. (2018) are expected to hamper the excitonic properties, while the extension of partially self-consistent schemes to the time domain is not straightforward. Therefore, benchmarks of different methods in and out of equilibrium will yield valuable insights.

In this work, we study a two-band semiconductor model in one dimension (the extended Falicov-Kimball model), with virtual excitons induced by a local inter-band interaction. This simple model has all the ingredients needed for exploring exciton dynamics far from equilibrium, and highly accurate solutions can be obtained. In particular, we employ the infinite time-evolving block decimation (iTEBD) Vidal (2007) method, which – upon convergence – yields an essentially numerically exact description. Furthermore, we employ several methods within the NEGF framework, including time-dependent mean-field (tdMF) theory and the full treatment of the Kadanoff-Baym equations (KBEs).Stefanucci and Leeuwen (2013) The self-energy is treated in the second-Born approximation (2BA), which can capture polarization and exchange effects. Furthermore, we employ the generalized Kadanoff-Baym ansatz,Lipavsky et al. (1986) which reduces the computational cost considerably. While the iTEBD method is a numerically powerful and reliable method for one-dimensional systems, it is difficult to extend the method to more general setups such as higher dimensions and long-range interacting systems. In the present study, we use it to calculate benchmark results for the other methods. Such a systematic comparison for finite systems demonstrated the potential of the GKBA.Schlünzen et al. (2017) Here, we will show that the GKBA also performs well in extended systems.

Benchmarking these methods against each other, we systematically study the properties of excitons out of equilibrium and discuss the effects which require a treatment beyond the MF theory. In particular, we compare above-bandgap excitations to resonant excitations of the exciton. We show that in the latter case a moderately strong pulse can induce large coherent oscillations in the polarization, which can gradually decay and be regarded as a transient nonequilibrium excitonic phase. We systematically study the nature of collective modes in transient states after above-bandgap excitations and resonant excitations and find that in the latter case its nature is different from the normal exciton states in equilibrium. Combining the GKBA and the RPA-like approach with distributions obtained from GKBA, we reveal that the peculiar collective mode originates from the efficient creation of an inverted population at the edge of the valence and conduction band. We also study the cooling effects from the electron-phonon couplings and show the real time formation of the peculiar mode from the above-bandgap excitation and the build-up of an exciton condensation out-of equilibrium.

The paper is organized as follows. In Sec. II, we introduce our model and the methods (tdMF, 2BA, GKBA and iTEBD) used to study the time evolution of the model after photo-excitation. We also derive the expressions for the relevant susceptibilities. In Sec. III, we show the results of the simulations. Section III.1 presents the results in the linear response regime, while in Sec. III.2 we go beyond the linear response regime and discuss the difference between above-bandgap excitations and resonant excitations. In Sec. III.3, we consider the effects of cooling from the electron-phonon coupling. The conclusions of our study are summarized in Sec. IV.

II Formulation

II.1 Model

In this paper, we focus on a spinless two-band model,

[TABLE]

where the first term represents the kinetic energy

[TABLE]

Here $\langle i,j\rangle$ indicates a pair of nearest-neighbor sites, and $a=c,v$ indicates the orbitals. $c$ and $v$ stand for the conduction band and the valence band, respectively. $\hat{c}^{\dagger}$ is the electron creation operator, $J_{a}({\bf r}_{ij},t)$ the hopping parameter, ${\bf r}_{ij}$ is the spatial vector connecting site $j$ to site $i$, and $\Delta_{a}$ is the energy of orbital $a$. The electrons in the two bands interact via a local interaction

[TABLE]

where $\hat{n}_{i,a}=\hat{c}_{i,a}^{\dagger}\hat{c}_{i,a}$. The effect of an external field is partially included in $J_{a}({\bf r}_{ij},t)$ via the Peierls substitution

[TABLE]

where ${\bf A}(t)=-\int^{t}{\bf E}(\bar{t})d\bar{t}$ is the vector potential, ${\bf E}(t)$ is the electric field, and $q$ the charge of the electron. This term corresponds to the intraband acceleration. The third term is the dipole excitation, which represents the interband excitation,

[TABLE]

Here the dipole matrix $\mathbf{d}_{c,v}$ is local and $\hat{\bf P}$ is the dipole moment per site. We use the notation $\bar{a}=c$ ($\bar{a}=v$) for $a=v$ ($a=c$). In the following, we set the length of the primitive vector, $\hbar$ and $q$ to unity.

Assuming translational invariance, we define the operators in momentum space, $\hat{c}^{\dagger}_{\mathbf{k},a}=\frac{1}{\sqrt{N}}\sum_{\bf k}e^{i{\bf k}\cdot{\bf r}_{i}}\hat{c}^{\dagger}_{i,a}$. Here $N$ is the number of sites. With these operators, one can express the Hamiltonian as

[TABLE]

with

[TABLE]

Here $\epsilon_{a}({\bf k})=-\sum_{l}J_{a}({\bf r}_{l})e^{-i{\bf k}\cdot{\bf r}_{l}}$, where the sum runs over nearest-neighbor sites.

Next we introduce the single-particle density matrix as

[TABLE]

Note that $\rho_{{\rm loc},a,b}(t)\equiv\rho_{ia,ib}(t)=\frac{1}{N}\sum_{\bf k}\rho_{{\bf k},a,b}(t)$. We also use $\boldsymbol{\rho}_{\bf k}(t)$ to express the $2\times 2$ matrix with elements $\rho_{{\bf k},a,b}(t)$.

In the present study, we consider one-dimensional chains and assume that the dipole matrix is directed along the chain and that ${\bf d}_{c}^{*}={\bf d}_{v}$. The system is excited with Gaussian pulses with various excitation frequencies.

II.2 Methods

In order to study the nonequilibrium dynamics of this system, we use several different methods: tdMF, the 2BA, the GKBA implemented with the 2BA and the iTEBD. In the following, we briefly introduce these methods and discuss the corresponding susceptibilities.

In general, a linear function $\chi^{R}_{\rm BA}(t,t^{\prime})=-i\theta(t-t^{\prime})\langle[\hat{B}(t),\hat{A}(t^{\prime})]\rangle$ can be measured by exciting the system with a weak excitation, $\hat{H}_{\rm ex}=F_{\rm ex}(t)\hat{A}$ with $F_{\rm ex}(t)\propto\delta(t-t^{\prime})$, and observing the evolution of $\hat{B}$. This is how we measure linear functions in the following. If $\hat{A}=\sum_{ij}A_{ij}\hat{c}^{\dagger}_{j}\hat{c}_{i}$ and $\hat{B}=\sum_{ml}B_{ml}\hat{c}^{\dagger}_{l}\hat{c}_{m}$, the response function can be expressed as

[TABLE]

where $\chi^{R}_{ml,ij}(t,t^{\prime})$ is the retarded part of the function

[TABLE]

defined on the the Konstantinov-Perel’ contour $\mathcal{C}$,Konstantinov and Perel’ (1961); Stefanucci and Leeuwen (2013) which runs from time [math] to time $t_{\text{max}}$ along the real time axis, back to zero, and then to $-i\beta$ along the imaginary time (Matsubara) axis. ${\mathcal{T}}_{\mathcal{C}}$ is the contour ordering operator and $t,t^{\prime}\in\mathcal{C}$ refer to contour arguments.

In particular, we consider the response function for $\hat{A}=\rho_{\nu,j}\equiv\hat{\Psi}^{\dagger}_{j}\boldsymbol{\sigma}_{\nu}\hat{\Psi}_{j}$ and $\hat{B}=\rho_{\mu,i}\equiv\hat{\Psi}^{\dagger}_{i}\boldsymbol{\sigma}_{\mu}\hat{\Psi}_{i}$, which we denote by $\chi^{R}_{\mu\nu}(t-t^{\prime};{\bf r}_{ij})$ for a steady state. Here $\hat{\Psi}_{i}=[\hat{c}_{i,c}\;\;\hat{c}_{i,v}]^{T}$ and $\boldsymbol{\sigma}_{\mu}$ is a Pauli matrix. In momentum space this response function is expressed as $\chi^{R}_{\mu\nu}(\omega;{\bf q})=\sum_{l}\int dte^{i\omega t}\chi^{R}_{\mu\nu}(t;{\bf r}_{l})e^{-i{\bf q}\cdot{\bf r}_{l}}$. Here, $\chi^{R}_{11}$ corresponds to the polarization-polarization response function.

II.2.1 Time-dependent mean-field theory

In the tdMF theory, we consider the time evolution of the one-particle density matrix Eq. (8) under the MF Hamiltonian, which is self-consistently determined at each time. Assuming translational invariance, the MF Hamiltonian is

[TABLE]

with

[TABLE]

The time evolution of the density matrix follows from the van Neumann equation, $\partial_{t}\boldsymbol{\rho}_{\bf k}(t)=-i[{\bf h}_{{\rm MF},{\bf k}}(t),\boldsymbol{\rho}_{\bf k}(t)]$ and the MF effect is taken into account through $\boldsymbol{\rho}_{{\rm loc}}(t)=\frac{1}{N}\sum_{\bf k}\boldsymbol{\rho}_{{\bf k}}(t)$. We also note that the Hartree term shifts the positions of the bands after the excitation since the occupation in the two orbitals changes.

Now we consider the linear response of a steady solution in the MF dynamics assuming that the steady state does not break the symmetry of the Hamiltonian (the system remains in the normal state). Here a steady solution means a state which does not change under the MF time propagation. The equilibrium state is one example. The expression for $\chi^{R}_{\mu\nu}(\omega;{\bf q})$ evaluated by the direct time propagation within the tdMF is

[TABLE]

Here $\boldsymbol{\chi}$ indicates the $2\times 2$ matrix whose components are $\chi^{R}_{\mu\nu}$ with $\mu,\nu=1,2$, and $\boldsymbol{\Theta}={\rm diag}[-\frac{U}{2},-\frac{U}{2}]$. $\boldsymbol{\chi}_{0}$ is the response evaluated by the time evolution without updating the mean field, which can be expressed as

[TABLE]

Here $\boldsymbol{\mathcal{G}}_{\mathbf{k}}(t)$ is the MF Green’s function at the steady-state, which is expressed as

[TABLE]

with vanishing off-diagonal components, since we assume that the steady state is a normal state. $E_{a}({\bf k})$ is the energy of the electron in band $a$ with momentum ${\bf k}$ determined with the MF Hamiltonian, Eq. (12), for the density distribution $n_{a}({\bf k})$. The explicit expression of the Fourier transformation of $\boldsymbol{\chi}_{0}(t;{\bf q})$ is

[TABLE]

with $\boldsymbol{W}_{c}=\begin{bmatrix}1&0\\ 0&0\end{bmatrix}$ and $\boldsymbol{W}_{v}=\begin{bmatrix}0&0\\ 0&1\end{bmatrix}$. We note that by using the equilibrium distribution $n_{a}({\bf k})=(1+\exp(\beta E_{a}({\bf k})))^{-1}$, Eq. (13) reproduces the well-known RPA-type susceptibility in equilibrium, which consists of ladder diagrams, see Appendix A.

One can simplify Eq. (16) for ${\bf q}={\bf 0}$ by introducing

[TABLE]

This rotation makes the off-diagonal elements of $\boldsymbol{\gamma}$ and $\boldsymbol{\gamma}_{0}$ zero, while

[TABLE]

We note that for positive frequencies ($\omega>0$), $\gamma_{11}(\omega)$ and $\gamma_{0,11}(\omega)$ are featureless, while $\gamma_{22}(\omega)$ and $\gamma_{0,22}(\omega)$ are responsible for nontrivial features in $\boldsymbol{\chi}$ and $\boldsymbol{\chi}_{0}$. In particular, $\chi_{11}=\frac{1}{2}(\gamma_{22}+\gamma_{11})$ implies that $\chi_{11}$ and $\gamma_{22}$ exhibit similar structures.

II.2.2 Full Kadanoff-Baym formalism: Second-Born approximation

In order to investigate the out-of-equilibrium correlated dynamics beyond the tdMF approximation, higher-order scattering processes need to be taken into account. The NEGF framework provides a systematic and versatile approach for treating many-body effects in the time domain. Stefanucci and Leeuwen (2013); Aoki et al. (2014) We define the general Green’s function GF on the Konstantinov-Perel’ contour $\mathcal{C}$ as

[TABLE]

Adopting again the matrix notation, the GF obeys the equation of motion (Dyson equation)

[TABLE]

where $\delta_{\mathcal{C}}(t,t^{\prime})$ is a straightforward generalization of the Dirac delta function to the contour $\mathcal{C}$, while $\ast$ denotes the convolution along $\mathcal{C}$. Solving Eq. (20) is accomplished by projecting onto observable times by invoking the Langreth rules, yielding the KBEs.Stefanucci and Leeuwen (2013); Aoki et al. (2014) After solving the corresponding equilibrium state (Matsubara GF), the real-time evolution is governed by the KBEs. Since the MF self-energy $\boldsymbol{\Sigma}_{\rm HF}(t,t^{\prime})=\delta_{\mathcal{C}}(t,t^{\prime})({\bf h}_{\mathrm{Hartree}}(t)+{\bf h}_{\mathrm{Fock}}(t))$ is included in ${\bf h}_{\mathrm{MF}}(t)$, many-body effects beyond mean field are captured by the correlation self-energy $\boldsymbol{\Sigma}_{\mathrm{corr}}=\boldsymbol{\Sigma}_{\mathrm{corr}}[{\bf G}]$, which is a functional of the GF. In this work, we employ the 2BA, which corresponds to the second-order self-consistent weak-coupling approximation: $\boldsymbol{\Sigma}_{\mathrm{corr}}[{\bf G}]\approx\boldsymbol{\Sigma}^{\mathrm{2B}}[{\bf G}]$. The correlated parts of the self-energy consists of a direct and and an exchange part,

[TABLE]

For the interaction Hamiltonian (3), the direct contribution to the self-energy reads

[TABLE]

while the exchange part is given by

[TABLE]

While exchange effects captured by Eq. (II.2.2) vanish when the GFs do not have inter-orbital components, their impact onto the strongly driven dynamics is less clear. Therefore, we also compare results within the simplified 2BA (taking the direct contribution Eq. (II.2.2)) to the full 2BA. We denote the simplified 2BA as s2BA in the following.

Given the expression of the self-energy, one can evaluate the linear response functions by simulating the evolution after a weak delta-function field pulse. Using the real-space representation for convenience, the corresponding response function ($\chi_{ml,ij}(t,t^{\prime})$ in Eq. (10)) can be expressed as

[TABLE]

Here, ${\bf G}_{0}$ indicates the full Green’s function without the probe excitation, $F_{\rm ex}(t;i,j)$ is the strength of the external field proportional to $\hat{c}^{\dagger}_{j}\hat{c}_{i}$, $\frac{\delta_{\mathcal{C}}}{\delta_{\mathcal{C}}}$ is the functional derivative on the contour, and $\frac{\delta_{\mathcal{C}}\boldsymbol{\Sigma}[{\bf G}](\bar{t}_{1},\bar{t}_{2})}{\delta_{\mathcal{C}}F_{\rm ex}(t^{\prime};i,j)}$ the reducible vertex expressed as a functional derivative on the contour $\mathcal{C}$. The matrix ${\bf e}_{ij}$ is defined by $[{\bf e}_{ij}]_{kl}=\delta_{ik}\delta_{jl}$. The self-energy $\boldsymbol{\Sigma}[{\bf G}]$ entering Eq. (II.2.2) is the full self-energy $\boldsymbol{\Sigma}[{\bf G}]=\boldsymbol{\Sigma}_{\mathrm{HF}}[{\bf G}]+\boldsymbol{\Sigma}_{\mathrm{corr}}[{\bf G}]$. We note that the contribution from $\frac{\delta_{\mathcal{C}}\boldsymbol{\Sigma}^{F}[{\bf G}](\bar{t}_{1},\bar{t}_{2})}{\delta_{\mathcal{C}}F_{\rm ex}(t^{\prime};i,j)}$ leads to the ladder diagrams consisting of $G_{0}$. In other words, the response to the probe evaluated by only updating $\boldsymbol{\Sigma}^{F}[{\bf G}]$ in the Dyson equation and keeping $\boldsymbol{\Sigma}_{\rm corr}[{\bf G}]=\boldsymbol{\Sigma}_{\rm corr}[{\bf G}_{0}]$ corresponds to the ladder diagrams consisting of $G_{0}$. Hence, $\frac{\delta_{\mathcal{C}}\boldsymbol{\Sigma}_{\rm corr}[{\bf G}](\bar{t}_{1},\bar{t}_{2})}{\delta_{\mathcal{C}}F_{\rm ex}(t^{\prime};i,j)}$ generates diagrams beyond these ladder diagrams.

II.2.3 Generalized Kadanoff-Baym ansatz

The numerical cost of evaluating the full Kadanoff-Baym equations Eq. (20) scales as $\mathcal{O}(N\cdot N_{t}^{3})$, where $N_{t}$ is the number of time points used in the simulation, and it grows significantly for long propagation times. Employing the GKBA reduces the computational effort by one order of magnitude in $N_{t}$ and thus allows simulations up to considerably longer times. Furthermore, the GKBA has been shown to cure some deficiencies of the full KBE approach, especially for finite systems.Schlünzen et al. (2017) Systematic assessments in extended system are scarce,Schüler et al. (2019) which is one of the motivations for the present study.

Within the GKBA, the description is reduced to the time evolution of the single-particle density matrix. Given a self-consistent approximation to the self-energy ($\boldsymbol{\Sigma}=\boldsymbol{\Sigma}[{\bf G}]$), the equation of motion for the density matrix (transport equation) can be expressed as

[TABLE]

where the collision integral ${\bf I}^{<}_{\mathbf{k}}(t,t)$ is defined by

[TABLE]

Here, we consider the Keldysh contour, which starts from $t=-\infty$, in constrast to the Konstantinov-Perel’ contour used in the previous section. In the Keldysh formalism, correlations of the initial state $\boldsymbol{\rho}(t=0)$ are built in by adiabatic switching: at $t=-\infty$, the equilibrium density matrix is determined by the MF treatment, while correlation effects are gradually incorporated by replacing $\boldsymbol{\Sigma}_{\mathrm{corr},\mathbf{k}}(t,t^{\prime})\rightarrow f(t)f(t^{\prime})\boldsymbol{\Sigma}_{\mathrm{corr},\mathbf{k}}(t,t^{\prime})$ with a smooth switch-on function $f(t)$. However, Eqs. (25) and (II.2.3) are not closed in terms of $\boldsymbol{\rho}$ since, in principle, information on the whole two-time dependence of the GF enters the collision integral Eq. (II.2.3).

The idea of the GKBA is to approximate the Green’s functions (GF) in the collision integral by combining the information contained in the occupation (${\boldsymbol{\rho}}$) and the spectrum ($\tilde{{\bf G}}^{R},\tilde{{\bf G}}^{A}$) by introducing the following auxiliary GF:

[TABLE]

Here we determine $\tilde{{\bf G}}^{R/A}(t,t^{\prime})$ as the mean-field GF

[TABLE]

The GKBA attains a closed form for any choice of the self-energy upon replacing $\boldsymbol{\Sigma}[{\bf G}]\rightarrow\boldsymbol{\Sigma}[\tilde{{\bf G}}]$ and ${\bf G}\rightarrow\tilde{{\bf G}}$ in the collision integral Eq. (II.2.3). In the present paper, we use the full 2BA Eq. (21) as well as the simplified version which considers the direct contribution Eq. (II.2.2) only (s2BA).

We now roughly discuss the relation between the susceptibility evaluated by GKBA and the full KBE form as described in the previous section. As mentioned in the previous section, keeping $\boldsymbol{\Sigma}_{\rm corr}[{\bf G}]=\boldsymbol{\Sigma}_{\rm corr}[{\bf G}_{0}]$ in the full KBE corresponds to the ladder diagram in terms of the full GF ${\bf G}_{0}$, which is in contrast to the tdMF, whose ladder diagram consists of the MF GF. In the latter GF, the damping of quasi-particles is not included. In the language of the transport equation, Eq. (25), this corresponds to keeping $\boldsymbol{\Sigma}_{\rm corr}[{\bf G}]=\boldsymbol{\Sigma}_{\rm corr}[{\bf G}_{0}]$ but updating ${{\bf G}}$ in the collision integral. In the GKBA we approximately update ${\bf G}$ and $\boldsymbol{\Sigma}$ in the collision integral. Therefore, naively speaking, the corresponding susceptibility should include a) the effects of the ladder diagrams consisting of dressed Green’s function (more than the MF Green’s function) and b) the effects beyond the ladder diagrams.

II.2.4 iTEBD

In this subsection, we briefly explain the principle of iTEBD.Vidal (2007) This method can be applied for the time-dependent problems such as quench dynamics or laser driving in quantum spinBarmettler et al. (2009); Takayoshi et al. (2014, 2019) and fermion Bauer et al. (2015); Ono et al. (2016); Coulthard et al. (2017); Ono et al. (2017) systems. The advantage of iTEBD is that calculations without finite size effects are possible by assuming translational invariance of the system.

In one dimension, the quantum states can be represented as matrix product states (MPS). When the system has a translational symmetry, the MPS representation is also translationally invariant

[TABLE]

where $s_{i}$ represents the quantum state on the site $i$, and in the present system $s_{i}=0,1,2,3$ correspond to $(n_{iv},n_{ic})=(0,0),(1,0),(0,1),(1,1)$, respectively ($n_{iv},n_{ic}$ is the eigenvalue of $\hat{n}_{iv},\hat{n}_{ic}$). $\alpha_{i}$ is the suffix for the matrices, and the values in the diagonal matrix $\lambda_{\alpha_{i}}$ ($=\lambda_{\alpha_{i}\alpha_{i}}$) are singular values (also known as the entanglement spectrum) obtained from the Schmidt decomposition on the bond between the sites $i$ and $i+1$. The bipartition of the sites into $A$ and $B$ is for the purpose of the time evolution described below.

The initial state is $s_{i}=1$ for all $i$, and the MPS representation is given as $\lambda_{\alpha_{i}=1}^{AB(BA)}=1$ and $\Gamma_{\alpha_{i}=1,\alpha_{i+1}=1}^{A(B)}[s_{i}]=\delta_{s_{i}1}$, where the matrix dimension is 1. Next we write the Hamiltonian in the bipartite form as

[TABLE]

where

[TABLE]

Note that $\hat{H}^{AB(BA)}(t)$ only acts on the bond $AB(BA)$. Using the Trotter formula, the time evolution operator $\mathcal{U}(t,t+\Delta t)$ for an infinitesimal time interval from $t$ to $t+\Delta t$ is decomposed as

[TABLE]

We can consider $\mathcal{T}_{s_{i}s_{i+1};s_{i}^{\prime}s_{i+1}^{\prime}}\equiv e^{-i\hat{H}_{i}^{AB}(t+\frac{\Delta t}{2})\frac{\Delta t}{2}}$ as a two-site quantum gate, and the procedure of its application is as follows. We construct a large matrix

[TABLE]

and then perform the singular value decomposition,

[TABLE]

by regarding $(\alpha_{i-1},s_{i})$ and $(s_{i+1},\alpha_{i+1})$ as the row and column of the matrix, respectively. The number of updated singular values $\tilde{\lambda}_{\alpha_{i}^{\prime}}^{AB}$ is four times larger than that of $\lambda_{\alpha_{i}}^{AB}$ because $\tilde{\lambda}$ is obtained from the enlarged matrix $\Theta_{(\alpha_{i-1},s_{i});(s_{i+1},\alpha_{i+1})}$ ($s_{i}=0,1,2,3$). Since the dimension of the matrix increases by iterating the step, we fix a maximum dimension $M$ (called the truncation dimension) and only keep the $M$ largest singular values, truncating the rest when the matrix dimension exceeds $M$. The updated $\Gamma$ is constructed as

[TABLE]

The procedure is the same for the application of $\prod_{i\in B}e^{-i\hat{H}_{i}^{BA}(t+\frac{\Delta t}{2})\Delta t}$. By iterating the above update, we can calculate the time evolution of the system. The numerical error arises from the Trotter decomposition and the truncation, and the precision becomes better for larger $M$ and smaller $\Delta t$. In this paper, we set $M=1200$ and $\Delta t=0.01$ or 0.05 depending on the laser field. The expectation value of a single-site observable such as $\hat{n}_{i,v}$ and $\hat{P}_{i}$ (for the $A$ site) is calculated as

[TABLE]

where $*$ represents the complex conjugate. We also calculate the expectation value for the $B$ site in the same way and take the average of $A$ and $B$.

For the calculations of space-time correlation functions, we use TEBD for finite size systems instead of iTEBD because the application of the single site operator at the initial time $t_{0}$ breaks the spatially translational invariance. We prepare the $N$ ($=\mathrm{even}$) site system $r=-N/2+1,\ldots,N/2$, and apply the operator at the site $r=0$. The scheme for the time evolution of TEBD is the same as that of iTEBD. Hence the response functions are obtained directly

[TABLE]

where $\chi_{11}^{>}(t,t_{0};{\bf r})=-i\langle\hat{P}_{\bf r}(t)\hat{P}_{\bf 0}(t_{0})\rangle$ is the greater part of the contour function $\chi_{11}(t,t^{\prime};{\bf r})\equiv-i\langle\mathcal{T}_{\mathcal{C}}\hat{P}_{\bf r}(t)\hat{P}_{\bf 0}(t^{\prime})\rangle$. This quantity Eq. (29) reveals the structure of collective modes at finite momenta. $\chi_{11}^{>}(t,t_{0};{\bf r})$ can be calculated as follows. Since the initial state is $(n_{i,v},n_{i,c})=(1,0)$ for all $i$, the initial MPS is represented by one-dimensional matrix as stated above. For the equilibrium correlation function, we apply $\hat{P}_{\bf 0}$ to this state ($t_{0}=0$), and calculate the time evolution using the Hamiltonian without laser up to the time $t$. Then $\hat{P}_{\bf r}$ is applied and taking the inner product with the initial state (and the phase factor $e^{iE_{0}t}$, $E_{0}$ is the ground state energy). For the dynamical correlation function under the laser, we evolute from the initial MPS up to the time $t_{0}$ with the Hamiltonian with laser driving and obtain the state $|\Psi(t_{0})\rangle$. Then we evolute the two states $|\Psi(t_{0})\rangle$ and $\hat{P}_{\bf 0}|\Psi(t_{0})\rangle$ from $t_{0}$ to $t_{1}$ with the Hamiltonian under laser and apply $\hat{P}_{\bf r}$ only to the latter. $\chi_{11}^{>}(t,t_{0};{\bf r})$ is obtained as the inner product of these two states.

We note that $-\mathrm{Im}\underline{\chi}^{<}(\omega;{\bf q})$ in equilibrium at $T=0$ is exactly the same as $-\mathrm{Im}\chi^{R}(\omega;{\bf q})$ for $\omega>0$. In general, when the contribution from the lesser part of $\chi_{11}(t,t^{\prime};{\bf r})$ is small, $\underline{\chi}^{>}(\omega;{\bf q},t_{0})$ can be approximated with the Fourier component of the retarded part $\chi^{R}(\omega;{\bf q},t_{0})$. In practice we use a window function $F_{\rm window}(t;t_{0})$ in the Fourier transformation Eq. (29), $\chi_{11}^{>}\rightarrow\chi_{11}^{>}F_{\rm window}$, since TEBD can only access rather short times. Specifically, we use $F_{\rm window}(t;t_{0})=F_{\rm gauss}(t-t_{0},\sigma)$ with F_{\rm gauss}(t,\sigma)\equiv\exp\bigl{(}-\frac{t^{2}}{2\sigma^{2}}\bigl{)}.

III Results

In the following, we choose the hopping parameters as $J_{c}=1,J_{v}=-1$ and consider half-filling systems in the semiconductor regime (with a band gap $>0$). In this case the valence band is fully occupied in the ground state at $T=0$, which is our initial state. The single particle spectrum obtained by the MF theory becomes exact for this state, as discussed in Appendix. A. For the other parameters, we use $\Delta_{v}=-3.2$ and $\Delta_{c}=1.2$, and $U=2.0$, which corresponds to a direct gap semiconductor with the band gap ${\mathcal{E}}_{\rm gap}=2.4$ at $T=0$, see Fig. 1. The choice of these parameters is motivated by those of some TMDs, which are characterize by a binding energy of a few hundred meV and a gap energy of a few eV.Hong et al. (2014)

III.1 Linear response regime

We first discuss the excitons in the equilibrium system. The exciton state is a bound state of an electron in the conduction band and a hole in the valence band. When we denote the energy necessary to excite an exciton from the equilibrium state by ${\mathcal{E}}_{\rm ex}$, the exciton binding energy ${\mathcal{E}}_{b}$ can be expressed as ${\mathcal{E}}_{b}={\mathcal{E}}_{\rm gap}-{\mathcal{E}}_{\rm ex}$. To measure ${\mathcal{E}}_{\rm ex}$, we excite the system with a very weak and short pulse, which includes a wide range of frequency components, and measure the induced dynamics of the dipole moment $P$. The exciton energy ${\mathcal{E}}_{\rm ex}$ manifests itself as a well defined oscillation in this quantity, and thus can be obtained by the Fourier transformation of $P(t)$. In Fig. 2(a), we compare the ${\mathcal{E}}_{b}$ evaluated in the above way for different methods (s2BA, GKBA+s2BA,tdMF, iTEBD). The results match perfectly, since in the present case one can show that the MF dynamics (RPA-type response), the GKBA and 2B are exact, see Appendix A. (2BA and GKBA+2BA are also exact.) More specifically, the response function evaluated by Eq. (18) with the $T=0$ occupation becomes exact. In Fig. 2(b), we show the corresponding $\gamma_{0,22}(\omega)$. The imaginary part of $\gamma_{0,22}(\omega)$ is essentially zero below the band gap. (It is finite in the figure because we use $0^{+}=0.02$ in Eq. (18b) for the numerical evaluation.) The real part has a peak at ${\mathcal{E}}_{\rm gap}$, which is related to the imaginary part by the Kramers-Kronig relation. The crossing of $\gamma_{0,22}(\omega)$ and $-2/U$ at $\omega<{\mathcal{E}}_{\rm gap}$ leads to a peak structure in the imaginary part of $\gamma_{22}(\omega)$, which corresponds to the exciton. For the one dimensional case, one can analytically show that $-{\rm Re}\gamma_{22,0}(\omega)$ diverges $\propto\frac{1}{\sqrt{\omega-{\mathcal{E}}_{\rm gap}}}$ around ${\mathcal{E}}_{\rm gap}$ for $\omega<{\mathcal{E}}_{\rm gap}$ and that the binding energy ${\mathcal{E}}_{b}$ scales as $\frac{U^{2}}{4*(J_{c}-J_{v})}$ for small $U$. We also note that, as long as the ground state is semimetallic, the exciton binding energy is independent of ${\mathcal{E}}_{\rm gap}$ in the present case. One can see this from Eq. (18b). The change of the gap by $\Delta{\mathcal{E}}_{\rm gap}$ just shifts $\gamma_{0,22}(\omega)$ by $\Delta{\mathcal{E}}_{\rm gap}$. Hence the pole position in $\gamma_{22}(\omega)$ is also shifted by $\Delta{\mathcal{E}}_{\rm gap}$, and the binding energy does not change.

In Fig. 2(c), we show the spectrum of the linear response function $-{\rm Im}\chi_{11}(\omega;q)$ for $\Delta_{v}=-3.2,\Delta_{c}=1.2$ and $U=2.0$ at $T=0$ evaluated by the TEBD. One can see a dispersive band below the particle-hole continuum, which corresponds to the (virtual) exciton states and their dispersion.

III.2 Beyond linear response

Now, we discuss the time evolution of the system during and after a photo-excitation beyond the linear response regime. In the following, we use $\Delta_{v}=-3.2$ and $\Delta_{c}=1.2$, and $U=2.0$, which gives ${\mathcal{E}}_{\rm gap}^{\rm eq}=2.4$ and ${\mathcal{E}}_{\rm ex}^{\rm eq}=1.93$ in equilibrium at $T=0$. We apply the Gaussian pulse with

[TABLE]

Here F_{\rm gauss}(t,\sigma)(=\exp\bigl{(}-\frac{t^{2}}{2\sigma^{2}}\bigl{)}) is the envelope function and

[TABLE]

is a ramp-up function which ensures that the evolution of the field around $t=0$ is smooth. In the following, we use

[TABLE]

with $N_{\rm cyc}=10$ unless we mention the condition specifically. Here $N_{\rm cyc}$ is the number of cycles included within $[-3\sigma,3\sigma]$ of the Gaussian envelope. We will consider two cases, i) an excitation into the particle-hole continuum ($\Omega>{\mathcal{E}}_{\rm gap}^{\rm eq}$) and ii) a resonant excitation of the excitons ($\Omega={\mathcal{E}}_{\rm ex}^{\rm eq}$). The former case is depicted in Figs. 1 and  2(c) with green arrows, while the latter is shown with blue arrows. We note that in the case of strong excitations, ${\mathcal{E}}_{\rm ex}$ shifts away from its equilibrium value (${\mathcal{E}}_{\rm ex}^{\rm eq}$) during the pulse, so that for a fixed pulse frequency, the system eventually deviates from the resonant condition. With this excitation protocol, we are going to investigate how the exciton frequency ${\mathcal{E}}_{\rm ex}$, the binding energy ${\mathcal{E}}_{\rm b}$, and the single particle spectrum are affected by the photo-doping of the system.

In Fig. 3 we first show the GKBA+s2BA results for the time evolution of the number of electrons in the conduction band, the dipole moment, and the total energy after different excitations. For $\Omega=3.0>{\mathcal{E}}_{\rm gap}^{\rm eq}$ (left panels), the number of excited charge carriers increases with increasing field strength in this regime. In the absence of a field, the bands are decoupled and the Hamiltonian conserves the number of particles in the conduction and valence band, respectively, which is correctly captured by the GKBA. As expected, since $\Omega$ is far from the exciton frequency, there is no prominent oscillation observed after the pulse, which lasts up to $t\approx 20$. For $\Omega=1.9\simeq{\mathcal{E}}_{\rm ex}^{\rm eq}$ (right panels), one can observe a non-monotonic increase of $n_{c}$ as a function of time as well as the field strength. This can be understood as a Rabi oscillation between the ground state and the exciton state. After the excitation ($t\gtrsim 30$), one can observe strong coherent oscillations in $P(t)$ with some frequency $\omega_{\rm coh}$, which persist for a long time after the pulse. The damping speed of these oscillations is enhanced with increasing field strength. From the Fourier transformation of these oscillations, one finds $\omega_{\rm coh}=2.13$ ($\omega_{\rm coh}=2.12$) for $E_{0}=0.1$ ($E_{0}=0.2$) at $t=60$ (The frequency is a bit $\sim 0.04$ increased from just after the pulse.). These values exceed the exciton frequency in equilibrium and the renormalized gap energy ${\mathcal{E}}_{\rm gap}^{\text{ren}}=1.91$ (${\mathcal{E}}_{\rm gap}^{\text{ren}}=1.85$). Here ${\mathcal{E}}^{\rm ren}_{\rm gap}$ is extracted from the instantaneous MF hamiltonian ${\bf h}_{\rm HF}[\boldsymbol{\rho}](t)$. Note that when the amplitude of the polarization becomes small the contribution of the Fock term becomes negligible and ${\mathcal{E}}^{\rm ren}_{\rm gap}$ is mainly determined by the Hartree shift. (For smaller field amplitude $E_{0}$, the oscillation frequency is still smaller than the renormalized gap energy.) As demonstrated in Figs. 3(c) and 3(f), the total energy ($E_{\rm tot}$) is conserved after the excitation.

Now let us compare the results obtained by the different numerical methods. In Fig. 4, we compare the density of the conduction band electrons and the polarization obtained by s2BA, GKBA+s2BA, MF and iTEBD for $\Omega=1.9$. In all cases, the strong coherent oscillations in the polarization persist even after the pulse. Among the approximate methods (s2B, GKBA+s2B, MF), GKBA provides the results closest to those of iTEBD. The most important difference between the tdMF and the rest is the damping of the induced coherent oscillations. Although GKBA still underestimates the damping compared to iTEBD, we find that the estimation of the damping within the GKBA is quantitatively better for the stronger fields. The s2BA can also show the damping of oscillations but it is generally weaker compared to GKBA and for $E_{0}=0.1,0.2$ it is very weak, while 2BA and GKBA match better as we further increase the field strength. Importantly, the peculiar feature of the coherent oscillations induced by the resonant excitation can be observed in iTEBD. For example, within iTEBD $\omega_{\rm coh}$ is $2.05$, while ${\mathcal{E}}_{\rm gap}$ is $1.86$ for $E_{0}=0.1$ around $t=60$. (Since the direct evaluation of the single particle gap in nonequilibrium iTEBD calculations is difficult, we estimate ${\mathcal{E}}_{\rm gap}$ from the density of excited charges $n_{c}$ and the resulting Hartree shift.) We also compare 2BA, GKBA+2BA, s2BA and GKBA+s2BA in Appendix C, but, in the present setup, the exchange term does not result in a significant change in the evolution of $P$ nor systematically improve the results compared to s2BA and GKBA+s2BA. This comparison suggests that the GKBA captures well the relevant features of the dynamics of the extended systems and that it is a useful method for systematic studies due to its relatively cheap computational cost.

Let us now comment on the relation between the strong coherent oscillations observed here and results reported in previous works.Comte and Mahler (1986); Östreich and Schönhammer (1993); Perfetto et al. (2019, 2019); Eastham and Littlewood (2001); Szymańska et al. (2006) After the excitation, the Hamiltonian conserves the number of electrons and holes, respectively. If the excited charge carriers are cooled down due to some coupling to thermal baths, the steady state after the relaxation should be described by a thermal state of the original Hamiltonian (Eq. (1) without excitation) with two different effective chemical potentials for the conduction band ($\mu_{c}$) and valence band ($\mu_{v}$), $\hat{H}^{M}=\hat{H}(0)-\mu_{c}\hat{N}_{c}-\mu_{v}\hat{N}_{v}$.Perfetto et al. (2019); Eastham and Littlewood (2001); Szymańska et al. (2006) Here $\mu_{c}$ and $\mu_{v}$ are determined such that the number of electrons and holes is the same as that just after the excitation. Since $\hat{H}^{M}$ corresponds to the original Hamiltonian with a smaller band gap, it can exhibit an excitonic insulating (EI) phase (exciton condensation out of equilibrium).Perfetto et al. (2019); Szymańska et al. (2006) The time evolution of the system is however described by $\hat{H}$, not by $\hat{H}^{M}$, so that this state exhibits oscillations of the polarization (off-diagonal component of the density matrix) with frequency $|\mu_{c}-\mu_{v}|$, which is of the order of the band gap. It has recently been shown in an independent work (Ref. Perfetto et al., 2019) using the same model as considered here and tdMF, that such a state can be realized even without baths by using a suitable pulse shape. Thus, the strong coherent oscillations in the polarization observed here can be understood as a transient realization of a nonequilibrium EI phase and the frequency of the oscillations above the renormalized band gap can be attributed to the effective chemical potentials in the two bands. Still, we have to note that the state just after the excitation is not exactly the thermal equilibrium state of $\hat{H}^{M}$, and our beyond-MF simulation shows that the transient EI phase can decay because of the scattering between the excited carriers. In the following, we focus on the properties of the transient states characterized by a gradually decreasing polarization. The effects of a coupling to phonon baths, which results in the realization of an equilibrium state of $\hat{H}^{M}$, are discussed in Sec. III.3.

To study properties of the transient states, we perform a pump-probe simulation using GKBA + s2BA. Namely, in addition to the first strong pump field, we add a second weak pulse (probe pulse) with some time delay. The shape of the probe pulse is chosen as

[TABLE]

In the following we use $\sigma_{\rm probe}=0.5$ and $E_{\rm probe}=0.01$ and neglect the vector potential of the probe pulse. Then we measure the dipole moment $P(t)$ and calculate the difference between the results with and without a probe pulse at $t_{\rm probe}$,

[TABLE]

To identify frequencies of oscillations induced by the probe pulse at $t_{\rm probe}$, we perform a Fourier transformation with a window function,

[TABLE]

Here $F_{\rm window}(t;t_{\rm probe})=F_{\rm gauss}(t-t_{\rm probe},\sigma)$ and $\sigma=20.0$ is used in the following. This time dependent spectral function can reveal the excitation structure of the transient state around $t=t_{\rm probe}$, when the oscillations induced by the pump pulse is not large or slower than the characteristic frequency induced by the probe pulse. We call the peak in $\Delta P(\omega;t_{\rm probe})$ as $\omega^{*}_{\rm coh}$ in the following.

In Fig. 5 we show the results of these analyses for $\Omega=3.0$ and $\Omega=1.9$, respectively. For $\Omega=3.0$ (above band-gap excitation, left four panels), one finds that there is almost no change in $\Delta P(t;t_{\rm probe})$ and hence $\Delta P(\omega;t_{\rm probe})$ after the pump pulse. With increasing field strength, the oscillation frequency ($\omega^{*}_{\rm coh}$) becomes smaller and at the same time, the life time of the oscillation becomes shorter. After the excitation, the band gap is reduced because of the Hartree shift from the photo carriers. Still, the frequency of the oscillation is within the shifted band gap, and thus the situation is qualitatively similar to the exciton state in equilibrium. Thus we can regard $\omega^{*}_{\rm coh}$ as a renormalized exciton energy, $\mathcal{E}^{\rm ren}_{\rm ex}$. We note that within GKBA + s2BA, the renormalized binding energy, $\mathcal{E}_{\rm b}=\mathcal{E}^{\rm ren}_{\rm gap}-\mathcal{E}^{\rm ren}_{\rm ex}$, is slightly increased to $0.50$ ($0.55$) for $E_{0}=0.2$ ($E_{0}=0.3$) from the equilibrium value $\mathcal{E}_{\rm b,eq}=0.47$. However, within GKBA + 2BA, even though $\mathcal{E}^{\rm ren}_{\rm ex}(=\omega^{*}_{\rm coh})<\mathcal{E}^{\rm ren}_{\rm gap}$, $\mathcal{E}_{\rm b}=0.42$ for $E_{0}=0.2$. Whether the enhancement of $\mathcal{E}^{\rm ren}_{\rm b}$ is genuine or not is thus unclear. (iTEBD can only access short times for $\Omega=3.0$.)

For $\Omega=1.9$ (resonant excitation, right four panels), one observes a very different behavior from the case discussed above. Namely, the frequency of the induced oscillations ($\omega^{*}_{\rm coh}$) increases for small $E_{0}$ from ${\mathcal{E}}_{\rm ex}^{\rm eq}$ and decreases for larger $E_{0}$. More remarkably, the frequency can exceed the renormalized band gap unlike the normal exciton in equilibrium. We note that, when the nonequilibrium states induced by the pump pulse show strong oscillations, the signal induced by the probe pulse also follows these oscillations and $\omega^{*}_{\rm coh}$ becomes similar to $\omega_{\rm coh}$. Hence, the gradual shift of $\omega^{*}_{\rm coh}$ after the pulse for $E_{0}=0.2$ can be attributed to the shift of $\omega_{\rm coh}$ itself. When the amplitude of the oscillations induced by the pump pulse is damped and becomes small, $\omega_{\rm coh}$ and $\omega^{*}_{\rm coh}$ are essentially the same, since both oscillations can be regarded as a small perturbation around the state without the oscillations. As in the case of $\Omega=3.0$, the life-time of the oscillations becomes shorter with increasing field strength.

Since the exciton states should be strongly affected by the transient quasiparticle occupations, we study the time evolution of the momentum distribution of the charges ($n_{c}(k),n_{v}(k)$). Since $n_{c}(k)$ and $1-n_{v}(k)$ behave identically, we only show $n_{c}(k)$ in Fig. 6. For $\Omega=3.0$ [Figs. 6(a) and 6(b)], the charges are excited at finite momenta which correspond to $\Omega=E_{c}(k)-E_{v}(k)$. Even though there occurs a slight redistribution and the occupation around $k=0$ becomes nonzero, most of the excited charges remain at nonzero momentum, and after the pulse the nonthermal distribution function remains almost constant. This is qualitatively similar to the MF dynamics, even though the latter lacks scattering and the occupation around $k=0$ remains almost zero after the pulse, see Fig. 14 in Appendix B. The slow intra-band relaxation is a consequence of the one-dimensional setup we are using, which implies that the scattering between charges is strongly restricted because of the momentum conservation and the energy conservation. It is expected that if we use a higher-dimensional lattice or consider electron-phonon scattering, one can observe a faster thermalization/redistribution process. In Sec III.3 we will analyze the effects of electron-phonon couplings.

For $\Omega=1.9$ [Figs. 6(c) and 6(d)], the charges are directly excited around $k=0$. After the pulse, the distribution function remains almost unchanged. The comparison with the results from tdMF (Fig. 14 in Appendix B) shows that the redistribution of the population due to scattering is indeed captured by GKBA, which yields a smooth distribution as a function of momentum and leads to an increase of the occupation near $k=0$. In both simulations, a fast approach to a steady value is observed after the pump, which is consistent with a change of the oscillation frequency during or quickly after the excitation. We note that for $E_{0}=0.3$ the particles are broadly distributed in the momentum space compared to the case for $E_{0}=0.2$.

To understand the origin of the qualitatively different collective modes ($\omega^{*}_{\rm coh}$) in the transient state after the pump pulse, depending on the excitation frequency, we now perform an RPA-type analysis using the essentially steady momentum distribution after the pulse. The idea of this analysis is the following. First, we extract, the momentum distribution $n_{c}(k)$ and $n_{v}(k)$ after the pulse from the GKBA simulation. We then substitute these $n_{c}(k)$ and $n_{v}(k)$ (neglecting the interorbital components $\langle\hat{c}^{\dagger}_{c}\hat{c}_{v}\rangle,\langle\hat{c}^{\dagger}_{v}\hat{c}_{c}\rangle$) into Eqs. (16), (18b) and (13) to estimate the nonequilibrium susceptibility. We note that this approximation corresponds to the MF dynamics starting from the distribution given by $n_{c}(k)$ and $n_{v}(k)$ (without interorbital component), which is a steady-state solution of the MF equation of motion.

In Figs. 7 and  8, we show $\gamma_{0,22}(\omega)$ and $\gamma_{22}(\omega)$ for different pump frequencies and amplitudes. For $\Omega=3.0$ (Fig. 7), as we increase the field strength, more electrons are excited to the conduction band and the band gap becomes smaller because of the Hartree shift, see Eq. (12b). As a consequence, the edge of the imaginary part of $\gamma_{0,22}$ is shifted to lower energies and the peak at the edge is reduced because of the finite density of conduction electrons around $k=0$, see Eq. (18b). The electrons stuck at non-zero momentum appear in the imaginary part of $\gamma_{0,22}$ as a local minimum around $\omega=3.0$. Because the imaginary part of $\gamma_{0,22}$ is connected to the real part through the Kramers-Kronig relation, these features in the imaginary part lead to a shift of the peak and a reduction of the height of the peak in the real part. Still the peak in the real part in $\gamma_{0,22}$ is prominent, which leads to a well defined in-gap mode appearing in the imaginary part of $\gamma_{22}$, see Fig. 7(d). The exciton binding energy (the distance between the peak and the dashed line in Fig. 7(d)) is gradually reduced with increasing pulse amplitude, which reflects the reduction of the height of the peak in the real part of $\gamma_{0,22}$.

For $\Omega=1.9$ (Fig. 8), we observe a suppression of the band gap with increasing field strength. Different from the case of $\Omega=3.0$, the excited charges directly accumulate at the bottom of the conduction band around $k=0$. This produces a more drastic change in $\gamma_{0,22}$ and hence in $\gamma_{22}$. For $E_{0}=0.06$, the peak structure around the renormalized gap ${\mathcal{E}}^{\rm ren}_{\rm gap}$ is strongly suppressed in $\gamma_{0,22}(\omega)$, but there still exists a crossing between the real part of $\gamma_{0,22}(\omega)$ and $-2/U$, which leads to a well-defined in-gap state as in equilibrium, see Fig. 8(c,d). When we further increase $E_{0}$, a population inversion ($n_{c}(k)>n_{v}(k)$) occurs around $k=0$, which is reflected in the positive value of Im$\gamma_{0,22}(\omega)$ around ${\mathcal{E}}^{\rm ren}_{\rm gap}$. Because of this population inversion near $k=0$, the imaginary part of $\gamma_{0,22}(\omega)$ crosses zero at a certain energy, which we denote by $\omega^{*}$. This zero-crossing can lead to a peak in the real part of $\gamma_{22}(\omega)$. To show this let us approximate $\gamma_{0,22}(\omega)\simeq\alpha+i\beta(\omega-\omega^{*})$ around $\omega^{*}$. Using Eq. (18c),

[TABLE]

This expression features a pole at $\omega=\omega^{*}+i(1+\frac{U}{2}\alpha)/(\frac{U}{2}\beta)$. If $(1+\frac{U}{2}\alpha)/(\frac{U}{2}\beta)$ is small compared to $\omega^{*}$ and the range in which the linearization of $\gamma_{0,22}(\omega)$ is justified, one can see a clear peak in the real part of $\gamma_{22}(\omega)$ around $\omega^{*}$. This condition is indeed satisfied in the present case, see $E_{0}=0.08,0.2$ in Figs. 8(a)(b), where $\alpha\gtrsim-2/U$, so that we end up with a clear peak in the real part of $\gamma_{0}(\omega)$, see Figs. 8(c)(d).

Thus, the RPA-type analysis qualitatively reproduces the dependence of the frequency of the collective mode $\omega^{*}_{\rm coh}$ on the excitation condition. For the above-gap excitation with $\Omega=3.0$, $\omega^{*}_{\rm coh}$ stays smaller than $\mathcal{E}^{\rm ren}_{\rm gap}$, whose character is similar to that of excitons in equilibrium. On the other hand, the mode observed for $\Omega=1.9$ above the renormalized band gap is explained by the population inversion just at the bottom of band (large $\beta$) and the moderate excitation, which results in a minimum of the real part of $\gamma_{0,22}(\omega)$ close to $-2/U$. Since a resonant excitation at the equilibrium exciton energy can quickly induce such populations, its naturally result in the peculiar coherent mode with frequency $\omega^{*}$. Furthermore, the RPA-type analysis predicts that the appearance of a well-defined peak in the real part of the susceptibility $\boldsymbol{\chi}_{0}(\omega;{\bf q=0})$ instead of the imaginary part leads to a phase shift of the oscillation against the probe pulse.

We now directly check this change in the transient susceptibility within GKBA + s2BA. Using GKBA, we can estimate the transient susceptibility through the pump-probe simulation as

[TABLE]

with $\Delta P(\omega;t_{\rm probe})$ defined in Eq. (35) and $E_{\rm probe}(\omega;t_{\rm probe})=\int dte^{i\omega t}E_{\rm probe}(t;t_{\rm probe})$. We note that this corresponds to $\chi_{11}^{R}(\omega;{\bf q=0})$ in equilibrium when $E_{\rm probe}$ is very weak. Considering the fact that the system is oscillating, we calculate the average of $\chi_{\rm GKBA}(\omega;t_{\rm probe})$ over the time interval $145\leq t_{\rm probe}\leq 155$ ($\tilde{\chi}_{\rm GKBA}(\omega)$) and show the results in Fig. 9(a)(b).

For small $E_{0}$, there is a peak in the imaginary part of $\tilde{\chi}_{\rm GKBA}$, see $E_{0}=0.04$ as an example, while for large enough $E_{0}$, the peak appears in the real part of $\tilde{\chi}_{\rm GKBA}$, see e. g. $E_{0}=0.2$. This is consistent with the RPA-type analysis. In Fig. 9(c), we show the renormalized gap (evaluated at $t=150$) and the frequency of the induced oscillation $\omega^{*}_{\rm coh}$ evaluated by the peak position of $|\tilde{\chi}_{\rm GKBA}(\omega)|$ as a function of $E_{0}$. The relative magnitude of these quantities switches around $E_{0}=0.08$ but the weak- and strong-field regimes are smoothly connected (no singular behavior). In Fig. 9(c), we also show the phase of $-\tilde{\chi}_{\rm GKBA}(\omega)$ at $\omega=\omega^{*}_{\rm coh}$. Reflecting a peak in the imaginary part for small $E_{0}$ and the one in the real part for large $E_{0}$, the phase quickly changes from a value close to $1.5$ to one close to [math] near $E_{0}=0.08$.

Although the GKBA and the RPA-type analyses agree qualitatively, there are several differences between them. First, compared to GKBA, the RPA-type analysis shows a larger frequency of the collective mode and a more abrupt switching of the phase, Fig. 9(c). Second, GKBA predicts that the signal in the crossover region becomes larger and the peak becomes sharper compared to the result for larger values of $E_{0}$, which is opposite to the behavior found in the RPA-type analysis. We also note that for $E_{0}=0.17,0.19$, the RPA-type analysis predicts a positive weight at the peak in Re$\gamma_{22}$, which originates from the fact that $\gamma_{0,22}(\omega)$ becomes smaller than $-2/U$ at $\omega^{*}$. Hence, the phase of $-\tilde{\chi}_{\rm GKBA}(\omega^{*}_{\rm coh})$ takes a value near $-\pi$. In addition, the RPA-type analysis predicts an infinite life-time of the in-gap states, while in the GKBA analysis these states can decay.

These differences may be attributed to i) the absence of the effects of the interorbital components in the RPA-type analysis, ii) the fact that GKBA partially takes into account the finite life-time of the quasiparticles as well as the corrections beyond the ladder diagrams from the correlated part of the self-energy. Neglecting the effects of the off-diagonal part in the density matrix should not be justified when the induced oscillations are long-lived as in $E_{0}=0.06\sim 0.1$ at $\Omega=1.9$. Hence, the transition from the normal-exciton like oscillation to the peculiar collective mode above the band gap is not fully captured within the RPA-type analysis. As for ii), the finite lifetime of quasi-particles can lead to a decay of the excitons and hence a finite lifetime, while the vertex corrections for the response functions beyond the RPA-type diagrams can renormalize the frequency of the oscillations.

Finally, we show the momentum resolved correlation functions evaluated by TEBD, Eq. (29). In Fig. 10, we show $-\mathrm{Im}\underline{\chi}^{>}(\omega;{\bf q},t_{0})$ just after the resonant excitation ($\Omega=1.9$ and $E_{0}=0.1$), see Fig. 2(c) for the equilibrium result. In equilibrium, there is a single exciton band below the electron-hole continuum. The correlation function after the resonant excitation exhibits several well-defined bands. Around $q=0$, the sign of $-\mathrm{Im}\underline{\chi}^{>}(\omega;{\bf q},t_{0})$ changes around $\omega=2.05$, which is consistent with the behavior of $-{\rm Im}\tilde{\chi}^{R}_{\rm GKBA}$ when the peculiar mode is generated, see Fig. 9(b). Interestingly, the positive signal above $2.05$ evolves into a well-defined branch at finite momentum which has a different dispersion than the exciton branch.

III.3 Effects of electron-phonon coupling

So far we have studied the dynamics of pure electron systems. However, in practice, there are nonzero electron-phonon (el-ph) couplings and the excited charge carriers can be cooled down. Especially in semiconductors, the relaxation in the conduction band can occur on a few tenth to a few hundred of femtoseconds and thus plays an import role.Bar-Ad and Chemla (1997); Bányai et al. (1995) The efficiency of the cooling depends on the strength of the el-ph coupling and the phonon frequency. Here we study the cooling effects using the GKBA. Namely, in addition to the self-energy from the el-el interaction, Eq. (II.2.2), we add the self-energy representing the el-ph coupling at the level of the Migdal approximation:

[TABLE]

Here, $D_{0}(t,t^{\prime})$ denotes the phonon GF and we assumed that the phonons are locally coupled to the densities of each band on each site. Since the total density per site is fixed, with this type of coupling, no dynamics of the phonon displacement is induced after the excitation, and the Hartree-like (Ehrenfest) term can be ignored. Moreover, such a coupling to the phonon bath does not change the symmetry of the Hamiltonian so that the number of excited charges is conserved after the pulse. We fix the phonon propagator to the equilibrium one (no feedback to the phonon subsystems), such that the phonons act as a heat bath. The phonon GF is obtained by Fourier transforming $D^{\gtrless}_{0}(t,t^{\prime})=\int d\omega/(2\pi)D^{\gtrless}_{0}(\omega)e^{-i\omega(t-t^{\prime})}$ and the fluctuation-dissipation theorem $D^{>}_{0}(\omega)=-iN_{\mathrm{B}}(\omega)B(\omega)$, $D^{<}_{0}(\omega)=-i[N_{\mathrm{B}}(\omega)+1]B(\omega)$ ($N_{\mathrm{B}}(\omega)$ is the Bose distribution). Here we consider the Ohmic spectrum

[TABLE]

with cutoff frequency $\omega_{c}$.

In Fig. 11, we show the evolution of $n_{k,c}$ and the results of the pump-probe simulation for the excitation above the gap ($\Omega=3.0$) with different excitation strength. One can see the relaxation of the excited carriers from finite momentum toward $k=0$, which was absent in the case without electron-phonon coupling. Reflecting the time evolution of the momentum distribution, the frequency of the collective mode induced by the probe field gradually increases. For the weaker excitation, the frequency of the collective mode remains below the band gap, while, for sufficiently strong excitations, at some point in time the frequency exceeds the renormalized band gap. The latter result is very similar to the resonant excitation case without el-ph coupling, where the photo electrons (holes) are directly created at the bottom (top) of the conduction (valence) band. The present calculation shows that, with the cooling induced by the el-ph coupling and for sufficiently strong excitation, the peculiar collective mode can also be induced by above band-gap excitations.

As discussed in Sec. III.2, when the system is completely relaxed after the excitation, the steady state reached should be described by the original Hamiltonian (Eq. (1) without excitation) with two different chemical potentials for the conduction band ($\mu_{c}$) and valence band ($\mu_{v}$), $\hat{H}^{M}=\hat{H}(0)-\mu_{c}\hat{N}_{c}-\mu_{v}\hat{N}_{v}$. Perfetto et al. (2019); Eastham and Littlewood (2001) Since the ground state of such a Hamiltonian can be an excitonic insulating (EI) phase, one can expect the appearance of large amplitude persistent oscillations of the polarization (nonequilibrium EI phase). In Sec. III.2, we showed that the resonant excitation can create a transient state close to such an EI phase, consistent with the recent results in Ref. Perfetto et al., 2019. Here we show that, in the presence of el-ph coupling, after the initial decay of the transient EI state, the nonequilibrium EI phase is recovered due to the cooling effect. As a result, large-amplitude persistent oscillations of the polarization reappear at long times.

In Fig. 12, we compare the time evolution with and without the phonon bath for the strong resonant excitation, which generates the excited electrons near the $\Gamma$ point. For the present field strength, the polarization damps quickly after the pulse in both cases. However, in the presence of the phonon bath, the polarization recovers after some time and exhibits persistent oscillations, which suggests the exciton condensation induced by the cooling of the excited charges. In Fig. 12(c), we show the evolution of the kinetic and total energies. The phonon bath gradually reduces the total energy. After the pulse, the kinetic energy also gradually decreases, but it starts to increase when the signal of the coherence of polarization starts to recover. This is consistent with the Bardeen-Cooper-Schrieffer (BCS) scenario, since the ordered state lowers the interaction energy at the cost of increasing the kinetic energy. In Fig. 12(d), we show the evolution of the momentum distribution of the conduction band electrons. One can clearly see that the electrons are more concentrated around the $\Gamma$ point compared to Fig. 6(d). Slow oscillations in the density distributions set in around $t=120$, where the polarization starts to be enhanced. These oscillations become less prominent in later times, which suggests the system approaches a steady state.

Finally, let us comment on a few points. First, we expect the emergence of the exciton condensation even in the case of the off-resonant excitation $\Omega=3.0$ if we simulate up to long enough times. Since the charges are excited to higher energy, it requires several scatterings with phonons for them to relax to the Gamma point. Second, a similar cooling and resultant condensation of excitons is expected for other types of el-ph couplings as long as the coupling does not break the symmetry. However, if the coupling is not to the total density on a given site, one cannot ignore the phonon displacement and the resulting change of the electron energy levels due to the Ehrenfest term. Since the phonon displacement is expected to damp and approach some steady value, the steady state of the electrons will be determined by $\hat{H}^{M}$, taking into account the change of the energy levels due to the phonon displacement in a self-consistent manner.

IV Summary and Conclusion

We have studied the fate of excitons in photo-excited semi-conductors using a spinless two band model in one dimension and different numerical methods; the tdMF, the 2BA, the GKBA implemented with 2BA and the iTEBD method. In the linear response regime at $T=0$, all these methods produce the exact linear response functions. Hence the exciton energies (${\mathcal{E}}_{\rm ex}$) can be accurately measured from the long-lived oscillations in the dipole moment after a weak excitation. Beyond the linear response regime, strong coherent oscillations in the polarization can be induced by resonant excitations. In particular, peculiar coherent oscillations characterized by a frequency larger than the semiconductor gap emerges for a properly chosen excitation strength. This behavior was also confirmed by the iTEBD simulation. We pointed out that these oscillations can be understood as a signature of the transient emergence of a nonequilibrium excitonic condensate,Östreich and Schönhammer (1993); Perfetto et al. (2019) which can be gradually suppressed as time evolves by interaction effects beyond mean field. Although 2BA and GKBA show such a suppression they underestimate it compared to the iTEBD reference data. Still, GKBA captures relevant properties of the coherent oscillations and provides the best agreement with iTEBD among these approximate methods.

Focusing mainly on the GKBA results, we have closely analyzed the collective modes ($\omega^{*}_{\rm coh}$) of the transient stats after resonant and above-band-gap excitations using a numerical pump-probe simulation. In the latter case, $\mathcal{E}_{\text{ex}}$(=$\omega^{*}_{\rm coh}$) is reduced mainly because of the photo-induced Hartree shift, but the exciton binding energy remains positive and thus the situation in the photo-doped state is qualitatively similar to an equilibrium state with reduced gap. For resonant excitations, $\omega^{*}_{\rm coh}$ tends to be increased from the equilibrium value of exciton $\mathcal{E}^{\rm eq}_{\text{ex}}$. When the excitation is weak, $\mathcal{E}^{\rm ren}_{\text{gap}}-\omega^{*}_{\rm coh}$ is still positive. On the other hand, for stronger excitations, it can become negative but the mode induced by the probe pulse is still well-defined. We revealed the origin of this characteristic behavior using the RPA-type expression of the susceptibility and the nonequilibirum distributions from the GKBA analysis. In particular, the peculiar mode characterized by the frequency above the band gap ($\omega^{*}_{\rm coh}>\mathcal{E}^{\rm ren}_{\text{gap}}$) originates from the photo-induced population inversion accompanied by a moderate number of excited charges and the sharp accumulation of electrons (holes) at the edge of the conduction (valence) band. The energy of this mode is determined by the energy up to which the photo-doped band is populated. We also studied the cooling effect from the electron-phonon coupling within GKBA. Because of the cooling of excited carriers, the frequency of the collective mode evolves in time. We demonstrated that the efficient cooling of excited carriers and a sufficient amount of photo-doping can induce the peculiar mode above the band gap even after above-gap excitations. We also simulated the build-up of the nonequilibrium exciton condensation in the phonon-cooled photo-doped state.

In the present study, we focused on a simplified model to benchmark the reliability of the methods and to explore potentially interesting phenomena. Our study shows that GKBA essentially captures the relevant physics, which enables systematic analyses for extend systems at a reasonable computational cost. In the future, it would be important and interesting to study the time evolution of excitons and charge distributions using more realistic models within GKBA. In addition, GKBA may be also be useful to study the real-time dynamics associated with the condensation of excitons or exciton polaritons out of equilibrium. The condensation problem has so far been mainly addressed with steady-state formalisms. A more realistic model study would provide microscopic and detailed insights into the various nonequilibrium phenomena observed in transition metal chalcogenides as well as semiconductors in cavities.

Acknowledgements.

YM, MS and PW were supported by ERC Consolidator Grant No. 724103 and the Swiss National Science Foundation via NCCR Marvel. The calculations were run on the Beo05 cluster at the University of Fribourg. M. S. thanks the Alexander von Humboldt Foundation for its support with a Feodor Lynen scholarship.

Appendix A RPA-type analysis in the linear response regime

For completeness, we provide a proof that the MF dynamics in the linear response regime is exact in the present model at $T=0$. To this end, we consider the linear response in terms of the nonequilibrium Green’s function (GF) formalism.Aoki et al. (2014); Stefanucci and Leeuwen (2013) The electron GF is defined on the Konstantinov-Perel’ contour (${\mathcal{C}}$) as in Eq. (19). We also introduce the correlation function on the contour as

[TABLE]

where $\mu,\nu=0,1$ and $\hat{\rho}_{\mu,{\bf q}}=\frac{1}{\sqrt{N}}\sum_{i}e^{-i{\bf q}\cdot{\bf r_{i}}}\hat{\rho}_{\mu,i}$ with $\rho_{\mu,i}\equiv\hat{\Psi}^{\dagger}_{i}\boldsymbol{\sigma}_{\mu}\hat{\Psi}_{i}$ and $\hat{\Psi}_{i}=[\hat{c}_{i,c}\;\;\hat{c}_{i,v}]^{T}$. The retarded part of this function is the susceptibility (the response function). At $T=0$ in the present model, the state with the valence band fully occupied ($\equiv|\Phi_{0}\rangle$) is the ground state when the band gap is sufficiently large. Therefore,

[TABLE]

Here $t\prec t^{\prime}$ indicates that $t^{\prime}$ appears later than $t$ in terms of the contour ordering. In addition, the single particle Green’s function within the MF theory is exact at $T=0$ in this model, since $\hat{c}^{\dagger}_{c,\bf k}|\Phi_{0}\rangle$ and $\hat{c}_{v,\bf k}|\Phi_{0}\rangle$ are also eigenstates and the corresponding eigen energies (measured from the ground state energy) are $\epsilon_{c}({\bf k})+U+D_{c}$ and $-\epsilon_{v}({\bf k})-\Delta_{v}$.

Now we consider the diagrammatic expression for $\chi_{\mu\nu}(t,t^{\prime};{\bf q})$ in terms of the full electron Green’s functions. The expression consists of a) the ladder diagrams (Fig. 13(a)), b) diagrams which include ring diagrams of the type shown in Fig. 13(b), and c) the ladder-like diagrams, which include at least one crossing of the interaction lines (Fig. 13(c)). However, one can show that the contributions from b) and c) are zero at $T=0$ in the present model, because of Eq. (41). A ring diagram consists of either $G_{cc}(t,t^{\prime})$ or $G_{vv}(t,t^{\prime})$, since $G_{cv}(t,t^{\prime})=G_{vc}(t,t^{\prime})=0$. When we write the time of the vertices on the ring as $t_{1}(=t_{N+1}),t_{2},..t_{N}$, both $t_{i+1}\succ t_{i}$ and $t_{i+1}\prec t_{i}$ must appear because of the periodic boundary condition. Hence, according to Eq. (41), the ring contribution should always vanish. For the ladder-like diagrams, let us write the times of the vertices in the lower lines as $t_{1},t_{2},,,t_{N}$ and those on the upper lines as $t^{\prime}_{1},t^{\prime}_{2},,,t^{\prime}_{N}$. The elements of $[t_{1},t_{2},,,t_{N}]$ and those of $[t^{\prime}_{1},t^{\prime}_{2},,,t^{\prime}_{N}]$ are identical since the interaction is instantaneous. To get a nonzero value for the lower part one needs $t\succ t_{1}\succ t_{2}\succ\ldots\succ t_{N}\succ t^{\prime}$ or $t\prec t_{1}\prec t_{2}\prec,,,\prec t_{N}\prec t^{\prime}$, while a nonzero upper part requires $t\succ t^{\prime}_{1}\succ t^{\prime}_{2}\succ,,,\succ t^{\prime}_{N}\succ t^{\prime}$ or $t\prec t^{\prime}_{1}\prec t^{\prime}_{2}\prec\ldots\prec t^{\prime}_{N}\prec t^{\prime}$. In a ladder with crossed interaction lines these two conditions cannot be simultaneously satisfied, so that the contributions from diagrams of the type shown in Fig. 13(c) also vanish. Therefore, only ladder diagrams can give a nonzero contribution to $\chi_{\mu\nu}$.

One can show that the summation of all the ladder diagrams leads to Eq. (13) with Eq. (14), where $\mathcal{G}_{k}$ is the exact equilibrium Green’s function at $T=0$. Hence, the susceptibility evaluated from the MF dynamics at $T=0$ is exact.

Alternatively, one can use the expression of the response function following Eq. (II.2.2)Stefanucci and Leeuwen (2013); Murakami et al. (2016a, b) to show that the response function obtained by the tdMF, s2BA and 2BA is exact. The ladder diagrams originate from the functional derivative of the Fock diagram ($\frac{\delta_{\mathcal{C}}\Sigma_{F}[G]}{\delta_{\mathcal{C}}F_{\rm ex}}$). Since the Fock term Eq. 12c is included in all of these methods, the corresponding susceptibility also includes the ladder diagrams. On the other hand, the equilibrium Green’s functions described by these approximations are exact. (From Eq. (41), the correlated part of the self-energy should be zero at $T=0$.) Hence the diagrams other than the ladder diagrams in the susceptibility vanish for the same reason as discussed above. Therefore, tdMF, s2BA and 2BA also produce the exact response functions at $T=0$ in this model.

As for the GKBA, $\rho_{cc,{\bf k}}=0,\rho_{vv,{\bf k}}=1,\rho_{cv,{\bf k}}=0,\rho_{cv,{\bf k}}=0$ yields $\boldsymbol{\Sigma}^{</>}_{\rm corr}[\tilde{{\bf G}}]=0$. Hence this is a steady state solution and the adiabatic switching of the interaction leads to this state. When the excitation with the off-diagonal field ($\sum_{i,a}\hat{c}^{\dagger}_{i,a}\hat{c}_{i,\bar{a}}$) is applied to this ground state, the linear response in $\hat{c}^{\dagger}_{i,b}\hat{c}_{i,b}$ should be zero, since the corresponding linear response function ($\simeq\langle\hat{c}^{\dagger}_{i,b}(t)\hat{c}_{i,b}(t)\hat{c}^{\dagger}_{i,a}(0)\hat{c}_{i,\bar{a}}(0)\rangle$) is zero due to the conservation of particles in each orbital. Hence the response of $\rho_{k,cc}$ and $\rho_{k,vv}$ against the field starts from $\mathcal{O}(E_{0}^{2})$. Therefore, $\tilde{G}_{cc}^{<},\tilde{G}_{vv}^{>}=\mathcal{O}(E_{0}^{2})$ and $\tilde{G}_{cv},\tilde{G}_{vc}=\mathcal{O}(E_{0})$. Using these facts and directly evaluating $\boldsymbol{\Sigma}^{</>}_{\rm corr}[\tilde{{\bf G}}]$, one can show that all components in $\boldsymbol{\Sigma}_{\rm corr}$ behave as $\mathcal{O}(E_{0}^{2})$. Hence in the linear response regime, the collision integral is still zero and the time evolution is the same as in the MF theory.

Appendix B Momentum distribution from tdMF

In Fig. 14, we show the momentum distribution of the conduction-band electrons ($n_{c}(k)$) evaluated with tdMF. For the above-gap excitation ($\Omega=3.0$), the charges are excited at finite momentum and are stuck there after the excitation because of the absence of scattering in tdMF. For the resonant excitation, there emerges some finite occupation around $k=0$ and $k=\pm 1.0$. The occupation around $k=0$ corresponds to the direct excitation, while that around $k=\pm 1.0$ corresponds to absorption of two photons.

Appendix C Effects of the exchange term

In Fig. 15, we compare the time-evolutions described by s2BA, 2BA, GKBA+s2BA, GKBA + 2B and iTEBD for $\Omega=1.9$ to see the effect of the exchange term (Eq. (II.2.2)). While there are rather clear effects on the number of photo-carriers, the inclusion of the exchange term does not generally lead to a quantitative improvement of the results. As for the effect on the time evolution of the polarization, it seems less prominent and again there is no clear improvement associated with the exchange term.