Comparing the generalized Kadanoff-Baym ansatz with the full Kadanoff-Baym equations for an excitonic insulator out of equilibrium

We investigate out-of-equilibrium dynamics in an excitonic insulator (EI) with a finite momentum pairing perturbed by a laser-pulse excitation and a sudden coupling to fermionic baths. The transient dynamics of the excitonic order parameter is resolved using the full nonequilibrium Green's function approach and the generalized Kadanoff-Baym ansatz (GKBA) within the second-Born approximation. The comparison between the two approaches after a laser pulse excitation shows a good agreement in the weak and the intermediate photo-doping regime. In contrast, the laser-pulse dynamics resolved by the GKBA does not show a complete melting of the excitonic order after a strong excitation. Instead we observe persistent oscillations of the excitonic order parameter with a predominant frequency given by the renormalized equilibrium bandgap. This anomalous behavior can be overcome within the GKBA formalism by coupling to an external bath, which leads to a transition of the EI system towards the normal state. We analyze the long-time evolution of the system and distinguish decay timescales related to dephasing and thermalization.

Simulating these processes can be challenging since an accurate but computationally feasible theoretical description is required for simultaneously dealing with strong external fields, many-particle interactions, and transient effects. The nonequilibrium Green's function (NEGF) approach can address all these challenges [30][31][32]: It is not limited to weak driving or linear response only, the many-particle correlations can be systematically included by construction of self-energy diagrams, and the real-time Green's function gives access to time-dependent observables such as densities, currents, total energies, and spectral functions. The drawback is in the computational effort for solving the dynamical equations of motion for the Green's function, which scale with the number of timesteps cubed. A simplification to this issue was proposed already over 30 years ago in Ref. [33] by reducing the two-time-propagation of the Green's function to the time-propagation of a time-local density matrix via the generalized Kadanoff-Baym ansatz (GKBA), * riku.tuovinen@utu.fi thereby reducing the computational scaling to the number of timesteps squared. While this approach was acknowledged and used already in the 1990s [34][35][36][37][38], its recent revival [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56] has made it possible to combine the NEGF approach with ab initio descriptions of realistic atomic, molecular, and condensed matter systems [57][58][59][60][61]. Recent development has further allowed for an equivalent but more efficient representation of the GKBA time evolution with only a linear scaling in the number of timesteps [62][63][64].
In this work, we consider ultrafast many-particle correlations in an excitonic-insulator system acting as a prototypical ordered-phase material [1,7,10,65]. Out-ofequilibrium dynamics in such systems with a symmetrybroken ground state has been shown to be extremely sensitive to all the intricacies in the electronic and lattice structure [7,10,29,53,66]. We drive the system out of equilibrium in two ways: (1) by an external laser pulse, and (2) by coupling to fermionic baths. We compare the resolved dynamics for the NEGF between the full Kadanoff-Baym equations (KBE) and the computationally less expensive GKBA. We find that while the laser-pulse excitation introduces rich transient dynamics with predominant oscillations given by a renormalized bandgap, the GKBA description, in contrast to KBE, does not damp to a stationary solution. This can be attributed to narrow spectral features of the GKBA, the character of the approximation for the propagators, and correlation-induced damping in the KBE solution [67]. Coupling to fermionic baths instead opens up a natural decay channel for the GKBA description as well, and we observe clear damping and even a transition from the excitonic to the normal state. We further characterize the nature of this phase transition by identifying separate arXiv:2007.07801v1 [cond-mat.str-el] 15 Jul 2020 decay timescales.
The paper is organized as follows. In Section II we introduce the model system, and we outline the main equations of the NEGF and GKBA approach. The outof-equilibrium dynamics due to external laser pulses and coupling to fermionic baths are analyzed in Section III. In Section IV we summarize our conclusions and discuss future prospects.

A. Model for the excitonic insulator
We model the excitonic insulator (EI) by a two-band system of spinless fermions [7,53] kα are the annihilation (creation) operators for electrons with momentum k in band α ∈ {0, 1} labeling the two bands, and ∆ α is the associated crystal field leading to the bandgap ∆ ≡ |∆ 0 − ∆ 1 |. In practice, we consider a real-space structure of two one-dimensional chains with periodic boundary condition, see Fig. 1(a). Each of these two real-space structures result in each of the two bands as seen in Fig. 1(b). The creation and annihilation operators in momentum and real space are related byd mα , where mα labels the real-space lattice site m of the chain α ∈ {0, 1}. For one-dimensional chains with nearest-neighbor hopping J α the energy band dispersion is kα = 2J α cos k. In this picture, the crystal field ∆ α can be readily identified as the local on-site energy for the lattice points. In the real-space picture the Hamiltonian in Eq. (1) then readsĤ where the matrix elements are chosen such that for nearest neighbors in each chain h 0 mα,nα = J and for on-site h 0 mα,mα = ∆ α with ∆ 0(1) = − (+) ∆/2. For all calculations in the present work, we set J = −1 and calculate energies in units of |J| and times in units of |J| −1 .
The electron-electron interaction is taken as an interband Hubbard interaction [7,53] introducing a local density-density interaction of strength U for the electrons between the two bands. The electronelectron interaction is the origin of excitonic pairing between an electron pocket at k = 0 and a hole pocket at k = π, see Fig. 1(b). The excitonic insulator phase is determined by a finite order parameter d † (k+π)0d k1 = 0 which spontaneously breaks the conservation of charge in each of the bands and the spatial symmetry. We discuss the practical evaluation of the order parameter in Sec. II D. The pairing introduces a finite hybridization between the bands and opens a gap.
An external laser pulse driving the above system out of equilibrium is modeled by a direct transition between the two bands [7] where we set the pulse shape as a gaussian: ]e −4.6(t−tc) 2 /t 2 c of amplitude A, frequency ω, and centering t c = 2πn p /ω with n p being the number of optical cycles. Using the transformation introduced below Eq. (1) we rewrite also Eq. (4) in real space. Since the laser-pulse term couples the two bands at equal kpoints, using the property k∈[−π,π) e ik(m−n) = δ mn we obtain a straightforward replacement The total Hamiltonian for the above setup combining the kinetic, interaction, and external terms then readŝ From now on, we use matrix representations of these objects in terms of the one-particle states in the real-space basis {|mα }: mα|Ĥ 0 +Ĥ ext (t)|nβ = h mα,nβ (t) and mα|Ĥ int |nβ = v mα,nβ (t). While the interaction term itself is instantaneous, in Eq. (3), we allow the strength of it to be time-dependent to describe adiabatic switching, which we will discuss in Section III.

B. Time propagation of the nonequilibrium Green's functions
We employ the nonequilibrium Green's function (NEGF) method where the Kadanoff-Baym equations are propagated in time [30][31][32][68][69][70][71][72][73][74][75][76][77][78][79][80][81]. The key quantity in the formalism is the one-particle Green's function, which we write in the one-particle basis of our model system, where z, z are time coordinates on the Keldysh contour γ with the contour-time-ordering operator T γ . The contour γ has a forward branch, z = t − ∈ [0, ∞), and a backward branch, z = t + ∈ (∞, 0], on the real-time axis, and also a vertical branch, z = −iτ ∈ [0, −iβ] on the imaginary axis, see e.g. [31]. Here we set, without loss of generality, the contour starting point at zero on the real-time axis, z ≡ t = 0. The creation and annihilation operators are represented in the Heisenberg picture, and the ensemble average, denoted by · · · , is taken as a trace over the density matrix. The Green's function matrix G(z, z ) is the solution to the integro-differential equation of motion (in matrix form) where h(z) is the one-particle Hamiltonian for the system, δ(z, z ) is a delta function on the Keldysh contour, and Σ[G] is the self-energy kernel containing all the information about many-particle and embedding effects. The integration is performed over the Keldysh contour through the Langreth rules [82,83]. Depending on the contour-time arguments, (z, z ), the double-time functions appearing in Eq. (8) can be represented in components: lesser (<), greater (>), retarded (R), advanced (A), left ( ), right ( ), and Matsubara (M) [31]. The selfenergy kernel Σ[G] can be obtained from an underlying Φ-functional, Σ[G] = δΦ[G]/δG, to guarantee the satisfaction of various macroscopic conservation laws [84], provided that the equations of motion are solved selfconsistently [85][86][87][88].
The Green's function provides a direct access to system observables such as densities and currents of the outof-equilibrium system. In particular, we are interested in the time-dependent one-particle reduced density matrix (TD1RDM) given by the time-diagonal of the lesser Green's function, ρ(t) ≡ −iG < (t, t). At the equal-time limit on the real-time axis, z = t − , z → t + , we obtain from Eq. (8) and its adjoint [46,48,49,53] where we separated the self-energy, Σ(t, t ) ≡ Σ HF (t)δ(t, t ) + Σ corr (t, t ), in time-local Hartree-Fock (HF) and time-non-local correlation (corr) contributions, and we also introduced the collision integral in terms of the correlation part [46,48,49,53] We use the one-particle basis of our model system to write the self-energy at the HF level [53,79] (Σ HF ) mα,nβ (t) = δ mn δ αβ pζ v mα,pζ (t)ρ pζ,pζ (t) − v mα,nβ (t)ρ nβ,mα (t), (11) and the correlation self-energy at the second-order Born (2B) level [53,79] We note that since our model describes spinless fermions, the spin-degeneracy factor [32,79], typically written for the direct terms [first terms on the right-side of Eqs. (11) and (12)], is here simply 1.
The combination of the equation of motion in Eq. (8) and the expressions of the self-energies in Eqs. (11) and (12) represents a closed set of equations for the full solution based on KBE. We solve these equations using the numerical library NESSi [89]. In particular, we solve the problem in momentum space and use a suitable MPI parallelization over momentum points, see Ref. [7] for details. In the full KBE solution, the collision integral in Eq. (10) also includes the initial-correlation part on the imaginary branch of the Keldysh contour ∼ β 0 dτ Σ corr (t, τ )G (τ, t) [89]. From now on, we refer to this approach as 2B@KBE.
An alternative approach to close the equation of motion for ρ in Eq. (9) is to employ the GKBA approximation [33,41] where we denoted ρ < ≡ ρ and ρ > ≡ 1 − ρ, and we represent the retarded/advanced propagators at the HF level [41,44] with T being the chronological time-ordering operator. We then use Eq. (13) in Eqs. (12) and (10), and then solve for the TD1RDM in Eq. (9) by using a timestepping algorithm [53,88]. While the inclusion of initial correlations has been shown to be possible also within GKBA [52,54,90,91], here we adiabatically switch on the many-particle interactions and only include the collision integral in the form of Eq. (10). For efficient computation, we additionally use a recurrence relation for constructing Eq. (14) due to its group property [32,53], and we employ optimized matrix (tensor) operations for the construction of the 2B self-energy [92]. From now on, we refer to this approach as 2B@GKBA.

C. Inclusion of fermionic baths
So far we have considered isolated systems being exposed to external drives locally within the system. Now we add a contribution from a bath environment, e.g., a particle reservoir or a biased electrode, described by [31,93,94] where kλ labels the k-th basis function in the λ-th bath. The bath energy dispersion depends on the Keldysh contour time z [95,96] where µ is the equilibrium chemical potential and V λ (t) is a generic excitation, such as a bias voltage, taking place at z ≡ t = 0. The bath is coupled to the EI system by the coupling Hamiltonian [31,93,94] where J mα,kλ are the coupling matrix elements between the EI system and the bath, which in general also depend on the Keldysh contour time z. In this work, we consider the "partitioned approach" [97,98] where the systems are brought in contact at z ≡ t = 0. These contributions in Eqs. (15) and (17) are then added to the total Hamiltonian in Eq. (6). We consider electronic interactions only within the EI system. Hence, for a noninteracting bath environment the relevant Green's functions are given by [31,[93][94][95] where f (x) = 1/(e βx + 1) is the Fermi function at inverse temperature β, and we used f We may then readily write the retarded/advanced bath self-energy, which is completely specified by the bath and coupling Hamiltonians [31,[93][94][95], where we introduced ψ λ (t, t ) ≡ t t dtV λ (t) and [95] and we used the Cauchy relation for the relative-time Fourier transform of Eq. (18), 1/(ω−ε kλ ±iη) = P(1/(ω− ε kλ )) ∓ iπδ(ω − ε kλ ), with η being a positive infinitesimal and P denoting the principal value [31]. It is important to notice that the bath self-energy is represented in the basis of the EI system because it describes the effect of "embedding" the EI system into the bath environment. We now assume the frequency content of the bath selfenergy is much broader than the energy scales in the EI system, known as the wide-band approximation (WBA). This approximation is justified here as we are concentrating on very low-energy excitations within the EI system at which the bath density of states is practically featureless [99][100][101][102]. In the WBA, the level-width matrix becomes independent of frequency, Γ λ (ω) ≈ Γ λ , which means it becomes time-local. Then, also the real part of the self-energy in Eq. (21) vanishes due to Kramers-Kronig relations. Thus, the retarded/advanced bath selfenergy is obtained by further summing over the bath in- Similarly, we obtain for the lesser/greater bath selfenergy [93,95,103] Due to the WBA, the frequency integral in Eq. (24) as such is not convergent but we use a cutoff frequency, ω c , based on the physical band edge of the bath given by the bath energy dispersion: Since the retarded/advanced bath self-energy was obtained as a time-local contribution in Eq. (23), it can directly be included in the HF propagators in Eq. (14) [46,49] The lesser/greater component of the bath self-energy in Eq. (24), in contrast, appears in an additional collision integral [46,49] I bath (t) whose contribution is added to Eq. (10). Also, the GKBA of Eq. (13) is used for the lesser/greater Green's functions in Eq. (26).

D. Accessing physical observables
The TD1RDM, ρ(t), as a solution to Eq. (9) naturally contains the information about the single-particle density on its diagonal, but also time-dependent expectation values of any single-particle operatorÔ may be extracted using it by [104] In our model system, we consider excitonic pairing between an electron pocket of the upper band (around k = 0) and a hole pocket of the lower band (around k = ±π), see Fig. 1(b). In practice, this means that in the EI phase d † (k+π)0d k1 = 0. Therefore, we average this object over the reduced Brillouin zone (RBZ), k∈[−π/2,π/2) ≡ k , and define this as the excitonic order parameter [7,53] where N k is the number of k points in the RBZ, N is the total number of real-space lattice points, and we introduced where the limiting case applies for infinite lattice sites.
In practice, we evaluate the RBZ sum numerically, but in most cases already N = 20 corresponds to the sinc function fairly reasonably. On the second line of Eq. (28) we used the transformation of the field operators between momentum and real-space, which also results in the alternating sign, (−1) m = (e iπ ) m . Momentum-averaged band populations could be obtained similarly. The total energy in the system can be divided in three contributions: (1) single-particle (or kinetic) en- , where h includes the single-particle Hamiltonian and the external field; (2) HF energy E HF (t) = 1 2 Re Tr[Σ HF (t)ρ(t)] corresponding to the time-local part; and (3) correlation energy E corr (t) = − 1 2 Im Tr[I(t)] being the remaining part of the collision integral after removing the HF part [32]. While the effect of exchanging energy between the EI system and the external bath could be included in this description, we perform the energy considerations only for the isolated system. The total energy then reads We can further calculate energy absorption during some time interval by the difference where t final is, e.g., the total propagation time, and t initial the time when some external fields are being switched on. Alternatively, this could also be evaluated from a Hellmann-Feynman formula E abs = , since the field depends explicitly on time but the expectation value only implicitly.
The nonequilibrium spectral function is defined as [46] A which is a matrix in the one-particle states of our model system. It is important to note that the GKBA in Eq. (13) satisfies the exact condition G R − G A = G > − G < . We then calculate a spatio-temporal Fourier transformation of the nonequilibrium spectral function with respect to the real-space lattice coordinates and the relative-time coordinate τ ≡ t − t [104] where N is the total number of lattice points and T ≡ (t + t )/2 is the center-of-time coordinate. In practice, we evaluate it by setting T to half the total propagation time, when the relative-time coordinate τ spans the maximal range diagonally in the two-time plane. Eq.  full KBE approach. In the GKBA data we have used N = 24 as the total number of lattice points, hence the energy bands consist of discrete peaks, in contrast to the k-resolved KBE data in Fig. 2(d). In the limit of infinite number of lattice sites, these would produce the continuum energy-band structure of the EI system. In equilibrium we see the gap opening due to the excitonic condensate, see Fig. 2(b). The energy axis is adjusted with the equilibrium chemical potential to take the Hartree shift into account. We also see that the 2B@GKBA equilibrium spectral function, obtained via the adiabatic switching procedure to be discussed in the next Section, is very similar to the HF one: The density of states is modified slightly but the overall structure remains. Importantly, the 2B@KBE spectral features are more broadened compared to 2B@GKBA.

A. Driving with a laser pulse
For all calculations, we consider our system to be in the EI phase by setting ∆ = 1.4 and U = 3.5 [53]. In Fig. 3 we exemplify the generic procedure for the timedependent simulations. For the description of interactions at the HF level, the initial equilibrium state can be obtained with a separate time-independent calculation [53], and consequently the out-of-equilibrium behavior can readily be analyzed starting from t = 0. Here, we are mainly interested in the description of interactions at the 2B@GKBA level, going beyond the mean-field description. For this analysis, we first need to prepare the correlated equilibrium state. This can be obtained by an initial time evolution (t < t 0 ) without external fields but adiabatically switching on the many-particle interactions in the 2B@GKBA self-energies [53]. After this, the out-of-equilibrium behavior, due to a laser excitation for example, can be studied (t ≥ t 0 ). We note in passing that the preparation step may consume a considerable amount of computational time [53], and it would be highly attractive to apply some sort of a restart protocol, e.g. of Refs. [52,54,90,91], for a separate calculation starting at t = t 0 including the initially correlated state. However, we have experienced in numerous tests (not shown) for this procedure to result in non-stationary behavior. We suspect the EI system considered here, possessing a symmetry-broken ground state with nonzero coherences on the off-diagonals of the density matrix [53], may not provide an applicable equilibrium state, at least in the context of Ref. [52].
Let us first look at a concrete example of the time evolution at the HF or 2B@GKBA level. We fix the number of optical cycles in the laser pulse for all simulations n p = 2, cf. Eq. (4). In Fig. 3, we see that for the HF evolution the absolute value of the order parameter |φ(t)| remains constant without the applied field and it is substantially reduced and oscillating after the photo-excitation (A = 0.4, ω = 1.5). On the level of 2B@GKBA, the adiabatic switching procedure keeps the system in the EI phase, which is stationary without the applied field. This condition might change for different values of U and ∆ [53]. When we apply the the laser excitation the out-of-equilibrium dynamics is roughly similar in HF and 2B@GKBA: In 2B@GKBA the oscillation frequency is slightly increased compared to HF (see also Fig. 4(c) and the consequent discussion later on). Next, we will focus on the 2B@GKBA case and thoroughly analyze how the EI system's response depends on the laser excitation.
Stronger driving amplitude in the laser pulse expectedly makes the initial transient response stronger. This can be seen in Fig. 4(a) for t − t 0 6 |J| −1 where the excitonic order parameter decreases rapidly from its equilibrium value. This, however, does not mean the excitonic condensate would melt completely. Instead, the order parameter remains at an oscillatory but nonzero steady-state value after the laser pulse. The frequency of these steady-state oscillations is independent of the driving amplitude as can be seen from the Fourier spectra in Fig. 4(c) and corresponds to the amplitude mode excitations. The Fourier spectra are calculated using Blackman-window filtering [105]. As we increase the excitation strength, namely A ≥ 1, the order parameter after the photo-excitation is, somewhat counterintuitively, negligibly reduced. We will address this point more thoroughly later on.
The system expectedly responds more strongly to the resonant driving. This is seen in Fig. 4(b) where we find the system to be most in resonance with the driving frequency ω = 1.5. However, while the 2B@GKBA solution properly describes the resonance condition, it still retains its oscillatory character because of the lack of damping in the HF propagators. The oscillations of the excitonic order parameter after the laser excitation are independent of the laser frequency as can be seen from the Fourier spectra in Fig. 4(c). We also show the Fourier spectra of the HF data (cf. Fig. 3). As we saw already in Fig. 3 the oscillation frequency in 2B@GKBA is slightly increased compared to HF, from 2.8 to approximately 3. These values can be attributed to the equilibrium system parameter for the noninteracting bandgap ∆ = 1.4 as we see even harmonics with frequencies 2n∆ (with n a positive integer) in the HF spectrum. The oscillation can therefore be associated with the crystal field; even though the bandstructure gets modified due to the electron-electron interaction, cf. Fig. 2, the transient signatures include the remnants of the crystal field. We can verify this finding by breaking the symmetry of our lattice model by introducing a next-nearest-neighbor hopping J = J/2 (HF2 in Fig. 4(c)), in which case also the odd harmonics appear with frequencies (2n + 1)∆. In the 2B@GKBA data, the higher order harmonics are more suppressed while the basic resonant frequency, related to a renormalized equilibrium bandgap, remains clearly visible in all cases independent of the laser amplitude or frequency. We compare the 2B@GKBA solution to that of the full 2B@KBE in Fig. 5. In the weak excitation regime A 1, the excitonic order parameter is nonzero in the long time limit and its value roughly agrees between the 2B@GKBA and 2B@KBE results. However, the 2B@KBE solution shows a considerably stronger damping than the one of 2B@GKBA. This is due to the quasiparticle corrections beyond HF, in contrast to the form in Eq. (14), and the consequent correlation-induced damping [67]. For instance, if the driving frequency is slightly off-resonant, namely ω = 2.0, the narrow spectral window of 2B@GKBA does not capture as much of the weight as the more broadened 2B@KBE which damps towards a slightly different steady-state value. In case of the resonant driving ω = 1.5, the reduction of the order parameter is in an excellent agreement between the 2B@GKBA and 2B@KBE results. On the other hand, the dynamics is qualitatively different for strong excitation strengths A ≥ 1.0. While in the 2B@GKBA the order is negligibly reduced, it is completely melted for the 2B@KBE propagation scheme and the EI system undergoes a transition to the normal state consistent with the GW level description reported in Ref. [7].
The dependence on the driving amplitude presents the main difference between the 2B@GKBA and the 2B@KBE solution. Within the 2B@GKBA the steadystate value of the order parameter may depend nontrivially on the driving amplitude. For instance, for pulse frequency ω = 2 the order parameter is maximally reduced around A = 0.6 in Fig. 4(a). Higher amplitude pulses seem not to break the electron-hole pairs, keeping the excitonic order parameter roughly at its equilibrium value. This means that how the laser pulse get absorbed to the EI system depends strongly on the width of the spectral features, which are more narrow in 2B@GKBA than in 2B@KBE, see Fig. 2. We analyze this behavior more in detail in Fig. 6, where we show the energy absorption calculated using Eq. (31) as a function of the driving amplitude and frequency for both the 2B@GKBA and the full 2B@KBE solution. We have checked (not shown) that possible finite-size effects in 2B@GKBA are negligible as larger number of lattice sites in the EI model leads to qualitatively similar data. For both cases, we observe that for smaller driving amplitudes (A 1) the energy absorption is expectedly maximal around the resonant frequency ω = 1.5 related to the renormalized equilibrium bandgap, cf. Fig. 4. However, for 2B@GKBA, if we follow a line at fixed frequency, e.g., at ω = 1.5, we see that the energy absorption oscillates with the driving amplitude. This is not the case for the full 2B@KBE solution, where higher-amplitude pulses straightforwardly lead to larger absorption. For the 2B@KBE solution the moderately large electron-electron interaction U = 3.5 gives already considerable broadening, resulting in energy absorption and consequently melting of the excitonic condensate at any amplitude A 1.5 (cf. Fig. 4). On the other hand, we may conclude that the 2B@GKBA description is reasonable at weak fields close to resonance, but this picture breaks down at stronger fields off-resonance due to nonlinear absorption and higher order scattering mechanisms.
An interesting observation in the analysis of the absorbed energy is a softening of the absorption edge with an increased excitation strength, see the dashed regions in Fig. 6. For A 1.5 this onset of nonlinear absorption also seems consistent between 2B@GKBA and 2B@KBE. We can understand this phenomenon by analyzing a static problem with a constant dipolar matrix element. Because the form of the excitation in Eq. (4) introduces a direct dipolar transition matrix element, d † k1d k0 , it pushes the lowest and highest bands away from each other which, in turn, moves the backfolded bands in the middle closer to each other, cf. Fig. 2(a). The electron-electron interaction, on the other hand, introduces a further coupling between the bands in the middle, d † (k+π)0d k1 , leading to a competition between the excitonic order and the dipolar matrix element. We can verify this behavior by looking at the energy-and momentum-resolved spectral function in Fig. 7. In this calculation, we consider the equilibrium system supplemented with a constant dipolar transition A as in Eq. (4), which then shows how the bandstructure would be affected by this form of an excitation, in general. While these equilibrium spectral functions do not exactly correspond to the laser-pulse situation, it provides us with some insight on the underlying mechanism. We see the gap closing around A = 0.6, which is in this case the critical point where the equilibrium system transforms from the excitonic to the normal state. Higher transition amplitudes introduce simply a rigid shift of the bands away from each other when the electron-hole interaction is no longer binding them together. It would also be feasible to calculate the nonequilibrium spectral function due to the short laser-pulse excitation. However, due to the competing mechanisms and in contrast to Fig. 7, it would show a very rich and complex spectrum of multiple photon-assisted side bands, and as clear interpretation as in Fig. 7 would be challenging.

B. Coupling to fermionic baths
We now consider each lattice site of the two chains in our EI system to be coupled to two different baths with equal coupling strength J mα,kλ in Eq. (17). As the level width or tunneling rate Γ in Eq. (22) depends not only on the coupling strength but also on the bath energy dispersion, we investigate the role of bath coupling by directly varying the strength of Γ . The bath filling is modified by a bias V λ (t) in Eq. (16) which we set to a constant value − (+) V for the bath connected to the α = 0 (α = 1) chain of the EI system. For the bath environment we additionally fix β = 100 in Eq. (24). This effectively resembles a zero-temperature limit at which the adiabatic switching procedure is consistent.
The procedure for analyzing the dynamics induced by the bath coupling is similar to the laser pulse excitation in the previous subsection. We first prepare the correlated equilibrium state by the adiabatic switching procedure [53] and then suddenly bring the system in contact with the baths. The excitonic order parameter responds to this external perturbation as seen in Fig. 8. Also in this case, for a description of the electronic correlations at the HF level only, the bath coupling could be introduced without the preparation step, and the corresponding dynamics shows only a straightforward decay process depending on the coupling strength between the EI system and the baths, see Fig. 8. This decay behaviour is drastically modified when the electronic correlations are described at the 2B@GKBA level. Next, we will analyze this in detail by looking at the dynamics after the bath coupling at t 0 = 150 when varying (1) the bias, (2) the bath coupling duration, and (3) the bath coupling strength.
The bias changes the overall decay timescale of the excitonic condensate. In Fig. 9(a) we fix the bath coupling strength Γ = 0.1 and the bath coupling duration t bath = 10, and we show the excitonic order parameter dynamics when the bias is increased from V = ±0 to V = ±0.5U . The final state can have nonzero excitonic order if the energy injected by the bias is not large enough to break the electron-hole pairs completely. However, even the bath coupling itself without bias lowers the order parameter compared to the equilibrium value. The initial transient at t − t 0 < 2 |J| −1 is completely specified by the bath coupling strength, and the consequent decay dynamics depends on the bias. The bath coupling duration does not change the overall decay timescale of the excitonic condensate. In Fig. 9(b) we fix the bath coupling strength Γ = 0.1 and the bias V = ±0.2U , and we expose the EI system to the baths for varying durations. The initial transient on all the curves collapses onto one decay process described by the bath coupling strength and the bias, see also Fig. 9(a). The final state can also in this case have nonzero excitonic order if the bath exposure duration is short enough, but a transition from the EI state to a normal state is introduced for longer exposure durations.
Increasing the bath coupling strength, while keeping the bias and exposure duration fixed, makes the system undergo a faster decay process towards the normal state, see Fig. 9(c). This is understandable since stronger bath coupling directly influences the exponential decay timescale in Eq. (25). However, for weaker couplings the initial transient shows competing mechanisms for breaking and recombining electron-hole pairs. Interestingly, we also observe multiple exponential decay timescales which will analyze in detail next. We show the decay timescales of Fig. 9(c) separately in Fig. 10 in logarithmic scale, and we see clearly that the initial transient in all the cases is also here completely specified by the bath coupling strength. We thereby refer to this mechanism as dephasing [7,106,107]. A second exponential decay process can be seen when the bath coupling is strong enough to melt the excitonic condensate completely related to thermalization [7,106,107]. In this case, the bias was fixed to V = ±0, and the thermalization appears slower than dephasing. However, as we have seen in Fig. 9(a) the bias will affect the overall decay timescale, and increasing the bias can also make the thermalization faster than dephasing. We will look closer into this effect next.
In Fig. 11 we show the numerically extracted decay exponents from a wide selection of simulated decay processes with varying bias and coupling. We see that the dephasing timescale, τ de , remains roughly constant (given directly by the bath coupling strength 1/τ de ≈ Γ ) while the thermalization timescale, τ th , is affected by the bias. The trend here is consistent with Fig. 9(a) where we observed that higher bias results in faster decay. This is also similar to Ref. [7] where τ th reportedly grows with the excitation strength.
We can gain some more insight into these decay timescales by looking at the energy-and momentumresolved nonequilibrium spectral function in Fig. 12. Compared to the equilibrium spectral function in Fig. 2 the bath coupling expectedly modifies the spectral features drastically. In Fig. 12(a) we see that already with zero bias the coupled system's gap starts closing. For larger bias [ Fig. 12(b) and 12(c)] the system evidently transforms towards the normal state, cf. Fig. 2(a). It is also interesting to note that compared to the excitation in Fig. 7, the spectral properties in the case of bath coupling, Fig. 12, behave more straightforwardly as there seem to be no competing effects. This picture also translates into the clean decay dynamics of the excitonic condensate seen in Fig. 9 and the disentangled decay timescales seen in Figs. 10 and 11.

IV. CONCLUSION
We have considered the out-of-equilibrium dynamics in a prototypical ordered-phase material, namely the excitonic insulator. We have studied out-of-equilibrium conditions due to a laser-pulse excitation and coupling the EI system to a fermionic bath. The calculations based on the nonequilbrium Green's function and the generalized Kadanoff-Baym ansatz showed that the excited EI system may undergo a transition towards the normal state when coupled to a bath. However, the isolated EI system perturbed by a laser pulse showed persistent oscillations in the excitonic order parameter but the excitonic order was found to not melt completely. The analysis of the absorbed energy showed a good agreement between the GKBA and KBE in the weak photo-excitation regime. However, for strong excitations the GKBA underestimate the energy absorbed by the pulse.
The character of the dynamics of the EI system, whether excited by a laser pulse or coupled to a bath, was attributed to the narrow spectral features of the GKBA formalism where no proper thermalization channel was found to be present for isolated systems, at least on the level of Hartree-Fock propagators. The bath introduces a suitable decay channel, and we identified separate decay timescales for the excitonic order parameter related to dephasing and thermalization. While we have concentrated on the EI system, we expect our findings to also be general for other symmetry-broken or ordered-phase systems, including e.g., superconducting [9,11,12,108,109] or charge-density-wave order [110][111][112].
The present implementation of the interacting system embedded in a bath environment, and the subsequent solution of the dynamical equations of motion of the NEGF at the level of the GKBA allows for addressing simultaneously long timescales and large systems. For future work, we therefore highlight the possibility of investigating time-resolved quantum transport in relatively large junctions with electronic correlations [46,48,58]. In addition, addressing these effects could provide another route for strong indications of exciton condensation since enhanced tunneling currents in electron-hole double bilayer sheets of graphene and transition-metal dichalcogenide have recently been observed [113][114][115]. The GKBA approach for time-resolved quantum transport could also prove pivotal in, e.g., addressing transiently emerging topological phenomena in Majorana tunnel junctions [116] with longlasting characteristic current oscillations.