# Can a pure state remain pure in the Unruh effect?

Saptarshi Saha<EMAIL_ADDRESS>Department of Physical Sciences,
Indian Institute of Science Education and Research Kolkata, Mohanpur - 741
246, WB, India  Arpan Chatterjee<EMAIL_ADDRESS>Department of
Physical Sciences, Indian Institute of Science Education and Research Kolkata,
Mohanpur - 741 246, WB, India  Chiranjeeb Singha<EMAIL_ADDRESS>Chennai Mathematical Institute, H1, SIPCOT IT Park, Siruseri, Kelambakkam
603103, India

###### Abstract

A uniformly accelerated detector in an inertial vacuum undergoes an
unavoidable dissipation, and the final steady-state becomes thermal. However,
to attain such a mixed state, there is no bound for the acceleration of the
single atomic detector. Here we show that the scenario is entirely different
for two atoms with the same energy levels. There exists a critical limit of
the acceleration for two atomic detectors, below which the purity of a
particular initial state can be preserved. We observe that the generator of
the dissipative dynamics (Lindbladian) is invariant under a weak symmetry
transformation at this limit. Hence one of the eigenstates of the symmetry
operator is unchanged during the evolution. This kind of state is called a
quantum dark state, which is essentially a decoherence-free subspace. As a
consequence, the system becomes localized, and it can skip the Unruh
thermalization. Beyond the critical limit, the symmetry is explicitly broken.
Therefore our results suggest that the system goes through a first-order
dissipative phase transition from a localized to a thermal phase.

###### pacs:

04.62.+v, 04.60.Pp

_Introduction._ –The indistinguishability between the quantum vacuum
fluctuation and the thermal fluctuation is one of the fundamental discoveries
in the last century davies_scalar_1975 ; sciama_quantum_1981 . Several authors
introduced such a connection between the thermal bath and the quantum vacuum
as the thermalization theorem SEWELL_1980 ; SEWELL_1982 ; Kay_1985 . The
theorem tells that the inertial vacuum behaves as a thermal bath for a
uniformly accelerated observer. The temperature of the bath is proportional to
the acceleration of the detector. Later Unruh put this theorem in the context
of particle creation, which after him is known as the renowned Unruh effect
unruh_notes_1976 ; davies_scalar_1975 ; crispino_unruh_2008 ; BELL1983 . It is
also closely connected with the Hawking radiation in the near horizon region
hawking_black_1974 .
A uniformly accelerated two-level system (TLS) perceives the vacuum state of a
free scalar field as a thermal bath in which it undergoes a non-equilibrium
evolution. As a result, the system will reach a thermal state at equilibrium
sciama_quantum_1981 ; kolekar_indistinguishability_2015 . Therefore, one can
map this type of system to an open quantum system (OQS), where the effective
temperature of the effective thermal bath proportional to the acceleration
($\alpha$) of the TLS Benatti-unruh2004 ; Breuer_2002 ; Juarez2019 ; atom-
photon .
Thermalization is an irreversible process. Boltzmann first introduced the
concept of irreversibility using the classical collision model Kardar2007 .
The quantum analog of this problem is first described by Von Neumann
von_neumann_proof_2010 . Since then, the quantum thermalization problem is one
of the prestigious topics in modern physics nandkishore_many-body_2015 ;
abanin_colloquium_2019 . The thermal state can be identified as a mixed state
whose entropy is maximum. The initial memory is lost during this irreversible
evolution. No local conserved quantities survive at thermal equilibrium.
Hence, it is a non-integrable phase abanin_colloquium_2019 . Anderson first
investigated non-thermal features in an isolated random lattice due to
impurity. There exist several localized states, which skip the thermalization.
It is known as Anderson-localization anderson_absence_1958 ;
cutler_observation_1969 . These phases are also integrable as various local
conserved quantities persist at equilibrium abanin_colloquium_2019 . Many-body
localization, quantum scars are the other examples of the non-thermal phases
in the closed many quantum systems nandkishore_many-body_2015 . Also, in the
case of many-particle OQS, thermal to non-thermal phase transition occurs for
a critical value of a particular system parameter saha2021 . The steady-state
configuration suddenly changes at the critical point. It is further known as
the dissipative phase transition (DPT) in the literature manzano_symmetry_2014
; horstmann_noise-driven_2013 ; buca_note_2012 ; albert_symmetries_2014 ;
lieu_symmetry_2020 ; kessler_dissipative_2012 ; minganti_spectral_2018 .
Several experimental verifications support the existence of DPT, e.g., Dicke
quantum phase transition in optical cavity baumann_dicke_2010 .
As the temperature goes down, the quantum effect dominates the thermal
dissipation. So, at very low temperature, the spatially correlated thermal
bath behaves as a source of long-range quantum correlations, and the
cooperativity between the atoms arises kessler_dissipative_2012 ; saha2021 .
This behavior of the thermal bath is called the common-environment effect
benatti_environment_2003 . As a result, a persistent entanglement occurs in
the system of non-interacting spins in a thermal bath. So, the final steady-
state is a non-thermal state carmichael_analytical_1980 ; braun_creation_2002
. The whole dynamics is protected by several symmetry operators in the non-
thermal limit lieu_symmetry_2020 ; buca_note_2012 ; saha2021 . Therefore the
eigenstates of the symmetry operator have the possibility to remain unchanged
throughout the dynamics. It also predicts the preservation of the purity of
particular initial states. These dissipationless states are called dark states
buca_note_2012 ; saha2021 . Even for the markovian systems, they can carry the
initial memories in the equilibrium, which has huge implications in quantum
computation, and quantum optics fleischhauer_quantum_2002 . The steady-state
solution has a degeneracy at this limit. When the temperature increases, the
degeneracy is explicitly broken, and the cooperativity between the atoms is
lost, and the steady-state becomes thermal saha2021 .
Motivated from the above discussions of the temperature-driven phase
transition, we show that the acceleration-driven phase transition is possible
for a uniformly accelerated non-interacting two spin system in an inertial
vacuum. It is well known that, for a single atom ensemble, the final steady-
state is always thermal. Non-thermality arises when we have more than one atom
in the ensemble Benatti-unruh2004 . An inter-atomic distance-dependent
thermal-nonthermal signature of the Unruh effect in resonance Casimir-Polder
interaction (RCPI) was reported earlier Marino:2016rfa . Moreover, the energy
shift due to Casimir-Polder interaction (CPI) is very small compared to the
Zeeman levels, so it is very hard to observe experimentally Marino:2014rfa .
We observe that for such a system with fixed inter-atomic distance, the low
acceleration regime is entirely given a notion of non-thermality. The dynamics
are invariant under a weak symmetry transformation. The emergence of the two-
atomic dark state is the most important feature of this non-thermal phase, and
beyond the critical limit of the acceleration, the steady-state becomes
thermal. The system acts as a collection of single atom ensembles in this
phase. Although, this effect is completely different from the non-markovian
Unruh effect, where the non-thermality arises due to non-markovian effects
moustos_non-markovian_2017 ; sokolov_unruh_2020 . We also use the exact
experimental values estimated by Bell, and Leinass BELL1983  to give a proper
notion of the non-thermal regime, which can be proven by using the appropriate
instrumental setup in the future.
Throughout the paper, we follow the given notation. The four-vector is denoted
by ‘$x$,’ the time is represented by ‘$t$’ and the proper time by ‘$\tau$.’ We
denote the space co-ordinate by ‘${\bf{x}}$,’ $({\bf{x}}=X,Y,Z)$.

_Dynamics of the accelerated atoms through a massless scalar field._ – Here,
we consider a system of two spin half particles uniformly accelerated in
Minkowski spacetime and interacts with a free massless scalar field. The atoms
are following the hyperbolic trajectory with respect to the inertial observer.
The positions of the atoms in terms of proper time are written by,

$\displaystyle
t_{m}(\tau)=\frac{1}{\alpha}\sinh\alpha\tau,\,\,X_{m}=\frac{1}{\alpha}\cosh\alpha\tau,$
$\displaystyle
Y_{m}=0,\,\,Z_{a}(\tau)=Z_{a},\,\,Z_{b}(\tau)=Z_{b}.\,\,(m=\\{a,b\\}).$ (1)

We define $L$ to be the proper distance between the two atoms located at the
positions $(t_{a}(\tau),\,\bf{x}_{a}(\tau))$ and
$(t_{b}(\tau),\bf{x}_{b}(\tau))$. The Hamiltonian of the system+field is
expressed as Benatti-unruh2004 ,

$\displaystyle\mathcal{H}=\mathcal{H}_{\text{\tiny
S}}^{\circ}+\mathcal{H}_{\text{\tiny L}}^{\circ}+\mathcal{H}_{\text{\tiny
SL}}~{}.$ (2)

Throughout the paper, we use the natural units, $\hbar=c=1$. Here we assume
the two spins have the same Zeeman levels, and the free Hamiltonian of the
system is $\mathcal{H}_{\text{\tiny S}}^{\circ}=\mathcal{H}_{\text{\tiny
S}}^{\circ 1}+\mathcal{H}_{\text{\tiny S}}^{\circ
2}=\frac{1}{2}\omega_{\circ}\sigma_{3}^{(1)}+\frac{1}{2}\omega_{\circ}\sigma_{3}^{(2)}$
(The superscript in the Pauli matrices represent the atom number).
$\sigma_{i}$ are the Pauli matrices, and $\omega_{\circ}$ is the frequency of
Zeeman levels. The free Hamiltonian of the scalar field
$\mathcal{H}_{\text{\tiny L}}^{\circ}$ can be expressed as,
$\mathcal{H}_{\text{\tiny
L}}^{\circ}=\int\frac{d^{3}k}{(2\pi)^{3}}~{}\omega_{k}a^{\dagger}(k)a(k)~{}.$
Here, $\omega_{k}$ is the frequency of the free scalar field, and
$a^{\dagger},a$ are the creation and annihilation operator of the quantized
field, respectively. The coupling Hamiltonian between the atom-scalar field is
assumed to be Benatti-unruh2004 ,

$\displaystyle\mathcal{H}_{\text{\tiny
SL}}=\lambda\sum\limits_{\mu=0}^{3}\big{[}\sigma_{\mu}^{(1)}\otimes\phi_{\mu}(x_{1})+\sigma_{\mu}^{(2)}\otimes\phi_{\mu}(x_{2})\big{]}~{},$
(3)

where $\lambda$ is the coupling constant, $\phi$ represents the scalar field
and $x_{1},\,x_{2}$ are the individual trajectories of two atoms. Here,
$\phi_{\mu}(x)=\sum^{N}_{a=1}\big{[}\chi^{a}_{\mu}\phi^{-}(x)\,+\,(\chi^{a}_{\mu})^{\dagger}\phi^{+}(x)\big{]}~{}$.
$\phi^{\pm}(x)$ is the positive and negative field operator of the free scalar
field and $\chi^{a}_{\mu}$ are the corresponding complex coefficients. The
two-spin system and the scalar field are initially uncorrelated, so the total
initial density matrix can be written as
$\rho(0)=\rho_{s}(0)\otimes|0\rangle\langle 0|$. $|0\rangle$ is the vacuum
state of the scalar field and $\rho_{s}(0)$ is the initial density matrix of
the system. Under the weak-coupling approximation, using the _Von-Neumann-
Liouville_ equation in the interaction frame and taking the trace over the
field variable, we get the dynamical equation of the system, which is written
below Benatti-unruh2004 ,

$\displaystyle\frac{d\rho_{s}(\tau)}{d\tau}=-i\Big{[}\mathcal{H}_{lamb},\,\rho_{s}(\tau)\Big{]}+\mathcal{L}\big{(}\rho_{s}(\tau)\big{)}~{}.$
(4)

$\mathcal{L}\big{(}\rho_{s}(\tau)\big{)}$ is the dissipator of the master-
equation and $\mathcal{H}_{lamb}$ is known as the Lamb-shift Hamiltonian that
leads to the renormalization of Zeeman Hamiltonian. This Lamb shift gives rise
to the CPI between the spins Singha:2018vaj ; Chatterjee:2019kxg ;
Saha:2020rko ; Tian:2016uwp . The forms are given by,

$\displaystyle\mathcal{L}\big{(}\rho_{s}\big{)}$ $\displaystyle=$
$\displaystyle\sum\limits_{a,b=1}^{2}\sum\limits_{j,k=1}^{3}\gamma^{ab}_{jk}\Big{(}\sigma_{b}^{k}\rho_{s}\sigma_{a}^{j}-\frac{1}{2}\\{\sigma_{a}^{j}\sigma_{b}^{k},\rho_{s}\\}\Big{)}~{},$
(5) $\displaystyle\mathcal{H}_{lamb}$ $\displaystyle=$
$\displaystyle-\frac{i}{2}\sum\limits_{a,b=1}^{2}\sum\limits_{j,k=1}^{3}\mathcal{S}_{jk}^{ab}\sigma_{a}^{j}\sigma_{b}^{k}~{}.$
(6)

The field correlation time-scale is taken to be smaller than the relaxation
time-scale of the system atom-photon . The dynamics are completely positive
and trace preservation holds. Here, we assume $\chi^{a}_{\mu}$ satisfies
$\sum\limits^{N}_{a=1}\,\chi^{a}_{\mu}\big{(}\chi^{a}_{\nu}\big{)}^{\dagger}\,=\,\delta_{\mu\nu}.$
So the field correlation functions are diagonal i.e.
$G_{ij}^{ab}(x-y)=\delta_{ij}G^{ab}(x-y)$. The expressions of the field-
correlation function is obtained by,

$\displaystyle G^{ab}(\Delta\tau)$ $\displaystyle=$ $\displaystyle\langle
0|\hat{\Phi}(\tau,{\bf{x}}_{a})\hat{\Phi}(\tau^{\prime},{\bf{x}}_{b})|0\rangle~{}.$
(7)

Here $\Delta\tau=(\tau-\tau^{\prime})$. The correlation function is assumed to
be stationary. For a massless scalar field, in 4-dimensional Minkowski
spacetime, the standard expression of this function is written as
birrell_quantum_1982 ; chowdhury_unruh-dewitt_2019 ,

$\displaystyle G(x,x^{\prime})=-\frac{1}{4\pi^{2}R^{2}}~{},$ (8)

where
$R=\sqrt{(t-t^{\prime}-i\epsilon)^{2}-|{\bf{x}}-{\bf{x}}^{\prime}|^{2}}~{}$.
This is called the positive-frequency Wightman function. Here $i\epsilon$ is a
small constant in the complex plane and acts as a regulator to avoid the
divergence and make sure about the function’s analytic property over the
chosen complex plane. The Fourier transform along the trajectory of the atoms
is given below Hu2015 ,

$\displaystyle\mathcal{G}^{ab}(\omega_{\circ})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\frac{\omega_{\circ}}{1-e^{-2\pi\omega_{\circ}/\alpha}}f_{ab}(\alpha)~{}.$
(9)

Here,
$f_{ab}(\alpha)=\frac{\sin\Big{(}2\,\frac{\omega_{\circ}}{\alpha}\sinh^{-1}(\alpha\,L/2)\Big{)}}{L\,\omega_{\circ}\sqrt{1+L^{2}\,\alpha^{2}/4}}$
for $a\neq b$, and $f_{ab}(\alpha)=1$ for $a=b$. The explicit form of
_Kossakowski_ matrix, which is an essential component of the Lindbladian, is
obtained by Benatti-unruh2004 ,

$\displaystyle\gamma^{ab}_{jk}=A^{ab}\delta_{jk}-iB^{ab}\epsilon_{jkl}\delta_{3l}-A^{ab}\delta_{3k}\delta_{3l}~{}.$
(10)

$A^{ab}(B^{ab})$ are the even (odd) combination of the Fourier transform
$\mathcal{{G}}^{ab}$ for the positive and negative Zeeman frequency. For a
massless scalar field, the explicit forms of $A^{ab},B^{ab}$ are written as
Hu2015 ,

$\displaystyle
B^{ab}=\frac{\lambda^{2}\,\omega_{\circ}}{8\pi}f_{ab}(\alpha),\,A^{ab}=B^{ab}\times\coth(\pi\omega_{\circ}/\alpha).$
(11)

Here $a,b=\\{1,2\\}$. We use $f_{ab}(\alpha)$ as $f$ in the remaining part.

_Dissipative phase transition (DPT) and Unruh effect._ – The steady-state
solution of the QME gives the notion of the thermal equilibrium of the system.
In the Liouville space, the QME can be written as,
$\frac{d}{d\tau}\hat{\rho_{s}}=\hat{\mathcal{L}}\hat{\rho_{s}}$. Here,
$\hat{\rho_{s}}$ is a column vector of dimension $N^{2}\times 1$ and
$\hat{\mathcal{L}}$ is the Lindbladian superoperator of dimension $N^{2}\times
N^{2}$ ($N$ is the dimension of the Hilbert space). $\hat{\mathcal{L}}$ is a
non-Hermitian matrix albert_symmetries_2014 . The eigenvector corresponding to
the zero eigenvalues is the steady-state solution of the QME. DPT is defined
as, for a continuous changing of any system parameter, when the degeneracy of
the eigenstate of zero eigenvalues of $\mathcal{\hat{L}}$ is abruptly broken,
then the system makes a transition between two different phases. This
definition is quite similar to the quantum phase transition, where the ground
state degeneracy is broken at a critical limit Sachdev_1998 . In this case,
the steady-state solutions depend on the behavior of $f$. For a single
spin-$1/2$ system and $2\times 2$ matrix, the observables are
$(\sigma_{x},\sigma_{y},\sigma_{z})$ (one constraint-trace preservation) and
the corresponding equations are known as the _Bloch equations_ Bloch1953 .
Following the above argument for two spin systems, there exist fifteen
observables. Still, as we have symmetry in the system (the Zeeman levels of
the spins are the same), the number of observables is reduced to nine, which
can be written as saha2021 ,

Figure 1: (a) shows the numerical plot of $f(\alpha)$ vs. $\alpha$, where
$\omega_{\circ}=10^{14}$ Hz, $L=6\times 10^{-7}$ m. For the low value of
$\alpha$, $f=1$. Increasing the value of $\alpha$, we get $0\leq f<1$, and for
a large value of $\alpha$, $f=0$. At the transition point $\alpha\approx
2\times 10^{22}\,ms^{-2}$, the value of $f(\alpha)$ starts to fall from 1. (b)
shows the plot of the distribution of the eigenvalues of $\mathcal{\hat{L}}$
for two different values of $f$ (effect of lamb shift is neglected and
$A^{11},B^{11}$ is fixed, $A^{11}=4,B^{11}=1$). Here $p$ is the position of
the eigenvalues [1 to 16] in the increasing order. For $f=1$, there is a
degeneracy of the zero eigenvalues (denoted by the red “+” sign), as $f$
decreases ($f=0.8$), the degeneracy is broken (denoted by the blue “+” sign).
(c), (d), (e) show the plot of the first-order derivative of
$\\{M_{z},M_{zz},M_{c}\\}$ _w.r.t_ $\alpha$. At the transition point of
$\alpha=\alpha_{c}$, the function diverges. Therefore, for $f=1$, we get a
localized phase, and $f<1$, it is a thermal phase.

$\displaystyle M_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\text{Tr}_{s}\Big{(}[\sigma_{i}^{(1)}\otimes\mathbb{I}+\mathbb{I}\otimes\sigma_{i}^{(2)}]\rho_{s}\Big{)}~{},$
$\displaystyle M_{ii}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\text{Tr}_{s}\Big{(}[\sigma_{i}^{(1)}\otimes\sigma_{i}^{(2)}]\rho_{s}\Big{)}~{},$
$\displaystyle M_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\text{Tr}_{s}\Big{(}[\sigma_{i}^{(1)}\otimes\sigma_{j}^{(2)}+\sigma_{j}^{(1)}\otimes\sigma_{i}^{(2)}]\rho_{s}\Big{)}~{}.$
(12)

Although the choice is not unique. The Bloch-type equations for our case is
written as,

$\left[\begin{array}[]{c}\dot{M}_{z}\\\ \dot{M}_{zz}\\\
\dot{M}_{c}\end{array}\right]=\begin{bmatrix}-A^{11}&0&2B^{12}\\\
B^{11}/2&-2A^{11}&A^{12}\\\
-B^{12}/2&2A^{12}&-A^{11}\end{bmatrix}\left[\begin{array}[]{c}M_{z}\\\
M_{zz}\\\ M_{c}\end{array}\right]+\left[\begin{array}[]{c}B^{11}\\\ 0\\\
0\end{array}\right]~{}.$ (13)

Here $M_{c}=M_{xx}+M_{yy}$. Expectation value of the other observables goes to
zero at $\tau\to\infty$. Here, $f$ varies from $0\leq f\leq 1$. There exist a
critical value of $\alpha=\alpha_{c}$ for fixed $L$, when $f\to 1$. In this
limit, $A^{11}=A^{12}$, $B^{11}=B^{12}$. The Lindbladian super-operator
$\hat{\mathcal{L}}$ is invariant under the following unitary transformation.
$\hat{U}(\kappa)\hat{\mathcal{L}}\hat{U}^{\dagger}(\kappa)=\hat{\mathcal{L}}$
where $\kappa$ is real. We define, $\hat{U}(\kappa)=\exp{[-i\hat{D}\kappa]}$
and $\hat{D}=D\otimes\mathbb{I}-\mathbb{I}\otimes D^{T}$. Here,
$D=\sigma_{x}\otimes\sigma_{x}+\sigma_{y}\otimes\sigma_{y}+\sigma_{z}\otimes\sigma_{z}$.
$\hat{D}$ is precisely the weak symmetry super-operator, as it commutes with
the full $\hat{\mathcal{L}}$, unlike the strong symmetry, where the operator
commute with all the individual components of $\hat{\mathcal{L}}$
lieu_symmetry_2020 . As a consequences of the Noether theorem, there exists a
local conserved quantity, which is written as,

$\displaystyle\frac{d}{d\tau}(M_{xx}+M_{yy}+M_{zz})=0~{}.$ (14)

As we increase the number of spins, several two quantum local quantities will
be conserved (i.e.
$Tr_{s}(\\{\sum_{i}\sigma_{i}^{n}\otimes\sigma_{i}^{n+1})\rho_{s}\\}$). Here,
$n$ is the spin-number and $i=(x,y,z)$. It is called the localization
abanin_colloquium_2019 . Due to this symmetry the final steady-state has an
initial value dependence. In FIG. 1. (b), we numerically plot of the
distribution of the eigenvalues of $\mathcal{\hat{L}}$ in the increasing
order. $A^{11},B^{11}$ is chosen to be fixed. The contribution from the Lamb
shift is neglected, so, all values are either zero or real negative. For
$f=1$, the Lindbladian superoperator has two zero eigenvalue so the steady-
state is degenerate, and as $f$ decreases, the degeneracy is broken. The
steady-state solutions for $f=1$ is given by,

$\displaystyle M_{z}^{eq}$ $\displaystyle=$
$\displaystyle\frac{M_{\circ}(3+4(M_{2}+M_{3}))}{3+M_{\circ}^{2}}~{},$
$\displaystyle M_{c}^{eq}$ $\displaystyle=$
$\displaystyle-\frac{M_{\circ}^{2}-4(M_{2}+M_{3})}{2(3+M_{\circ}^{2})}~{}.$
(15)

Here, $M_{c}(\tau=0)=M_{2}$, $M_{zz}(\tau=0)=M_{3}$ and
$M_{\circ}=\tanh(\pi\omega_{\circ}/\alpha)$. The final steady-state is a non-
thermal one. This steady-state solutions are exactly matches with the previous
works by Benatti et al. Benatti-unruh2004 . The form of the symmetry operator
($D$) and the dipolar-coupling Hamiltonian are similar ficek_entangled_2002 .
Dicke first pointed out that the maximally entangled anti-symmetric state or
the singlet state [$|\psi\rangle=\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle)$]
which is one of the eigenstate of $D$, decays with a reduced spontaneous
emission rate in a common electro-magnetic field. This phenomena is known as
sub-radiance ficek_entangled_2002 ; Dicke1954 . Here also, if we choose the
value of $M_{2}=-1/4$ and $M_{3}=-1/2$, (which corresponds to the singlet
state), the final steady state solution is same as the initial
state,($M_{z}^{eq}=0$, $M_{zz}^{eq}=M_{2}$, $M_{c}^{eq}=M_{3}$). The purity of
such an initial state is conserved, $[\text{Tr}_{s}(\rho_{s}^{2})=1]$. Hence,
the dark states are also known as a sub-radiant state.
If $\alpha>\alpha_{c}$, then we get $0\leq f<1$. The Lindbladian is not
invariant under the above transformation. Hence the initial value-dependency
is gone and the steady-state is unique one. In this limit, the thermalization
theorem holds Takagi1986 . The steady-state solution is is given by,

$\displaystyle M_{z}^{eq}$ $\displaystyle=$ $\displaystyle M_{\circ}~{},$
$\displaystyle M_{zz}^{eq}$ $\displaystyle=$
$\displaystyle\frac{M_{\circ}^{2}}{4},~{}M_{c}^{eq}=0~{}.$ (16)

The thermal phase has no symmetries and no local conserved quantity. For a
relativistically uniform acceleration $\alpha$, the corresponding temperature
is given by, $kT=\hbar\alpha/2\pi c$. $k$ is Boltzmann constant, $\hbar$ is
Planck’s constant, and $T$ is temperature. Bell and Leinaas showed in a 3.1 km
circular accelerator in SPEAR at Stanford, and nearly 1200 k temperature can
be produced by circulating the ultra-relativistic electrons at a centripetal
acceleration of $\alpha=2.9\times 10^{23}m/s^{-2}$ BELL1983 . So we, vary
$\alpha$ from $(10^{21}\,ms^{-2}\to 10^{25}\,ms^{-2})$. Distance between the
atoms is chosen to be $L=6\times 10^{-7}m$, as the Lamb shift or the Casimir-
Polder interaction can be observed at a distance of $\mu m$ range Bressi2002
and Zeeman frequency is chosen as $\omega_{\circ}=10^{14}$ Hz. In principle,
localization in a thermal bath occurs at zero temperature. Several
experimental techniques like reservoir engineering, laser cooling ion-trap,
etc., are developed in the last few decades, which are used to enhance the
bath-correlation length at $T>0$ to trap the atoms Phillips1998 ; Poyatos1996
; Stellmer2013 ; Wineland1987 . Therefore, one can produce such a localized
phase even at non-zero temperatures. Similarly, the non-thermality will arise
at $\alpha=0$. Lowering the inter-atomic distance value and the increasing
value of Zeeman frequency will increase the field correlation length. As a
result, we can get a localized phase at $\alpha\neq 0$. These limits, we are
using here, are also experimentally justifiable BELL1983 .
For this choice of value, we numerically plot $f(\alpha)$ vs. $\alpha$. In
FIG. 1. (a), for $\alpha<2\times 10^{22}ms^{-2}$, we get $f=1$, which
signifies a localized phase and increasing $\alpha$, $f$ falls below 1, which
denotes the existence of thermal phase for higher acceleration. So we get two
different regions of $f$, which signifies two separate phases.
There exists an explicit symmetry breaking in the system. The first-order
derivative respect to $\alpha$, of the set of observables
$\mathcal{O}=(M_{z},M_{zz},M_{c})$ diverge at $\alpha\to\alpha_{c}$ (_i.e_ ,
$\lim_{\alpha\to\alpha_{c}}\frac{d\,\mathcal{O}}{d\alpha}=\infty~{}$). Beyond
the critical limit, the steady-state configuration abruptly changes.
Therefore, a first-order phase transition (from a localized to a thermal
phase) happens in the system. In FIG. 1. (c), (d), (e), we numerically show
that, the first-order derivative of $(M_{z},M_{zz},M_{c})$ diverges at
$\alpha_{c}\approx 2\times 10^{22}ms^{-2}$, which also predicts the existence
of two different phases. Non zero value of $M_{c}$ in the non-thermal phase
indicates the existence of field-induced entanglement in the system Benatti-
unruh2004 . Concurrence is a well-known measure of the entanglement between
the two spin systems Wootters1998 . The expression for the concurrence
($C(\rho_{s})$) is given by,
$C(\rho_{s})=\text{max}\\{0,\,4|M_{c}|-\sqrt{(1+4M_{zz})^{2}-4M_{z}^{2}}\\}$.
On the other hand, in the thermal phase, the entanglement vanishes
($M_{c}=0$), and the spins behave as a separate system.

_Conclusion._ – For a two-atom system accelerating through a massless scalar
field, the vacuum state of the field acts like a correlated bath. At the lower
values of acceleration ($f=1$, $\alpha<\alpha_{c}$), the system goes to a
localized phase, and the inertial vacuum acts as a common environment. There
exist a conserved quantity. Increase the value of acceleration, the symmetry
of the system is explicitly broken ($f<1$, $\alpha<\alpha_{c}$), the system
goes to a thermal phase. In this two limit ($f=1,f<1$), the steady-state
solution of the observables $(M_{z},M_{zz},M_{c})$ are different. Hence the
first-order derivative of the observables diverges at the transition point.
The Von-Neuman entropy $(-Tr_{s}(\rho_{s}\ln\rho_{s}))$ in the non-thermal
phase is not an extensive quantity saha2021 . In this phase, the dark states
do not evolve. Therefore, this kind of state is putting a limitation on the
thermalization theorem. When we increase the number of atoms in the ensemble,
the vacuum state from the point of view of the inertial atomic detector still
behaves as a canonical ensemble as a point of view of a uniformly accelerated
observer. But at the lower value of acceleration, the system is at a frozen
state, which is protected by several symmetry operators, so it cannot feel any
temperature. The principle of equal a priori probabilities do not hold here.
So, a pure dark state remains pure.

We wish to thank Dr. Rangeet Bhattacharyya and Pragna das for the fruitful
discussions. SS and AC acknowledge University Grant Commission (UGC) for
providing a research scholarship, and CS acknowledges Dr. Amitabh Virmani for
providing a visiting fellowship at Chennai Mathematical Institute.

## References

  * (1) P. C. W. Davies, Journal of Physics A: Mathematical and General 8, 609 (1975).
  * (2) D. W. Sciama, P. Candelas, and D. Deutsch, Advances in Physics 30, 327 (1981).
  * (3) G. L. Sewell , Phys. Lett. 79A, 23 (1980).
  * (4) G. L. Sewell , Ann. of Phys. 141, 201 (1982).
  * (5) B. S. Kay, Commun. Math. Phys. 100, 57 (1985).
  * (6) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * (7) L. C. B. Crispino, A. Higuchi, and G. E. A. Matsas, Rev. Mod. Phys. 80, 787 (2008).
  * (8) J. S. Bell and J. M. Leinaas, Nuclear Physics B 212, 131 (1983).
  * (9) S. W. Hawking, Nature 248, 5443 (1974).
  * (10) S. Kolekar and T. Padmanabhan, Class. Quantum Grav. 32, 202001 (2015).
  * (11) F. Benatti and R. Floreanini, Phys. Rev. A 70, 012112 (2004).
  * (12) H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, New York, 2002).
  * (13) B. A. Juárez-Aubry and D. Moustos, Phys. Rev. D 100, 025018 (2019).
  * (14) C. Cohen-Tannoudji, G. Grynberg, and J. Dupont-Roc, Atom-Photon Interactions: Basic Processes and Applications (Wiley, New York, 1992).
  * (15) M. Kardar, Statistical Physics of Particles (Cambridge University Press, Cambridge, 2007).
  * (16) J. von Neumann, EPJ H 35, 201 (2010).
  * (17) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * (18) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * (19) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * (20) M. Cutler and N. F. Mott, Phys. Rev. 181, 1336 (1969).
  * (21) S. Saha and R. Bhattacharyya, ArXiv:2104.08891 [Cond-Mat, Physics:Quant-Ph] (2021).
  * (22) D. Manzano and P. I. Hurtado, Phys. Rev. B 90, 125138 (2014).
  * (23) B. Horstmann, J. I. Cirac, and G. Giedke, Phys. Rev. A 87, 012108 (2013).
  * (24) B. Buča and T. Prosen, New J. Phys. 14, 073007 (2012).
  * (25) V. V. Albert and L. Jiang, Phys. Rev. A 89, 022118 (2014).
  * (26) S. Lieu, R. Belyansky, J. T. Young, R. Lundgren, V. V. Albert, and A. V. Gorshkov, Phys. Rev. Lett. 125, 240405 (2020).
  * (27) E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. Lukin, and J. I. Cirac, Phys. Rev. A 86, 012116 (2012).
  * (28) F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Phys. Rev. A 98, 042118 (2018).
  * (29) K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Nature 464, 7293 (2010).
  * (30) F. Benatti, R. Floreanini, and M. Piani, Phys. Rev. Lett. 91, 070402 (2003).
  * (31) H. J. Carmichael, J. Phys. B: Atom. Mol. Phys. 13, 3551 (1980).
  * (32) D. Braun, Phys. Rev. Lett. 89, 277901 (2002).
  * (33) D. Moustos and C. Anastopoulos, Phys. Rev. D 95, 025020 (2017).
  * (34) B. Sokolov, J. Louko, S. Maniscalco, and I. Vilja, Phys. Rev. D 101, 024047 (2020).
  * (35) M. Fleischhauer and M. D. Lukin, Phys. Rev. A 65, 022314 (2002).
  * (36) L. Rizzuto, M. Lattuca, J. Marino, A. Noto, S. Spagnolo, R. Passante, and W. Zhou, Phys. Rev. A 94, 012121 (2016).
  * (37) J. Marino, A. Noto, and R. Passante, Phys. Rev. Lett. 113, 020403 (2014).
  * (38) C. Singha, Mod. Phys. Lett. A 35, 1950356 (2019).
  * (39) A. Chatterjee, S. Saha, and C. Singha, EPL 130, 50004 (2020).
  * (40) S. Saha, C. Singha, and A. Chatterjee, Eur. Phys. J. C 81, 265 (2021).
  * (41) Z. Tian, J. Wang, J. Jing, and A. Dragan, Scientific Reports 6, 35222 (2016).
  * (42) N. D. Birrell and P. C. W. Davies, Quantum Fields in Curved Space, Reprint edition (Cambridge University Press, Cambridge, 1984).
  * (43) C. Chowdhury, A. Das, and B. R. Majhi, Eur. Phys. J. Plus 134, 65 (2019).
  * (44) J. Hu and H. Yu, Phys. Rev. A 91, 012327 (2015).
  * (45) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 2001).
  * (46) R. K. Wangsness and F. Bloch, Phys. Rev. 89, 728 (1953).
  * (47) Z. Ficek and R. Tanas, Physics Reports, 372, 369 (2002).
  * (48) R. H. Dicke, Phys. Rev. 93, 79 (1954).
  * (49) S. Takagi, Progress of Theoretical Physics Supplement 88, 1 (1986).
  * (50) G. Bressi, G. Carugno, R. Onofrio, and G. Ruoso, Phys. Rev. Lett. 88, 041804 (2002).
  * (51) W. D. Phillips, Rev. Mod. Phys. 70, 721 (1998).
  * (52) J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 77, 23 (1996).
  * (53) S. Stellmer, B. Pasquiou, R. Grimm and F. Schreck, Phys. Rev. Lett. 110, 26 (2013).
  * (54) D. J. Wineland, W. M. Itano, Physics Today, 40, 34 (1987).
  * (55) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).


# Cheers: a linear-scaling KBE+GKBA code

###### Abstract

The interaction of electrons with quantized phonons and photons underlies the
ultrafast dynamics of systems ranging from molecules to solids, giving rise to
a plethora of physical phenomena experimentally accessible using time-resolved
techniques. Green’s function methods offer an invaluable interpretation tool
since scattering mechanisms of growing complexity can be selectively
incorporated in the theory. cheers is a general-purpose nonequilibrium Green’s
function code that implements virtually all known many-body approximations and
is designed for first principles studies of ultrafast processes in molecular
and model solid state systems. The aims of generality, extensibility,
efficiency, and user friendliness of the code are achieved through the
underlying theory development and the use of modern software design practices.
Here, we motivate the necessity for the creation of such a code and overview
its design and capabilities.

###### keywords:

Nonequilibrium Green’s function theory, generalized Kadanoff-Baym Ansatz,
excited states

Y. Pavlyukh R. Tuovinen E. Perfetto G. Stefanucci

Y. Pavlyukh
Institute of Theoretical Physics, Faculty of Fundamental Problems of
Technology, Wroclaw University of Science and Technology, 50-370 Wroclaw,
Poland
<EMAIL_ADDRESS>
R. Tuovinen
Department of Physics, Nanoscience Center, P.O. Box 35, 40014 University of
Jyväskylä, Finland
E. Perfetto and G. Stefanucci
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca
Scientifica 1, 00133 Rome, Italy and INFN, Sezione di Roma Tor Vergata, Via
della Ricerca Scientifica 1, 00133 Rome, Italy

##  1 Introduction

Coherent electron dynamics in correlated materials is typically accompanied by
the interaction with bosonic particles and quasiparticles, such as phonons,
plasmons, charge density waves, and photons. Numerous proposals have been
pushed forward for creating materials with novel properties by putting them
into an entangled light-matter state [1]. Vastly different energy and length
scales, quantum aspects of the involved bosonic particles, and numerous
intertwined orders pose considerable challenges for theory. A scalable
theoretical method to model excitation and relaxation phenomena in correlated
many-body systems, reliable beyond the perturbative regime, is therefore
crucial to simulate and interpret experimental results and to suggest new
materials to study.

Many-body perturbation theory represents a systematic way to deal with inter-
particle correlations. In order to get access to the dynamical properties of
the system, the equations of motion (EOM) for the two-times electron and boson
Green’s functions (GF), hereafter referred to as the _nonequilibrium Green’s
function_ (NEGF) theory [2], must be propagated. The EOMs in this case are
known as the Kadanoff-Baym equations (KBE). The time non-locality of the
scattering term represents the major difficulty for the full two-times
propagation leading to the scaling that is at least cubic ($t_{\text{f}}^{3}$)
with the physical propagation time $t_{\text{f}}$ [3, 4, 5, 6, 7, 8, 9] making
it very difficult to resolve smaller energy scales associated, for instance,
with phonons [10]. Still, interesting progress has recently been made [11, 12,
13, 14]. The so-called generalized Kadanoff-Baym ansatz (GKBA) [15] allows one
to limit the propagation to time-diagonal, that is, to work with one-particle
density matrices rather than with two-times Green’s functions, preserving
concomitantly the description of inter-particle correlations. First-principles
implementations of the KBE+GKBA equations have been pioneered in atoms and
organic molecules in Refs. [16, 17, 18, 19, 20, 21, 22, 23] and 2D materials
[24, 25] using the cheers (a tool for the correlated hole-electron evolution
from real-time simulations) code [26].

KBE+GKBA became a competitive first-principles method in pure electronic case
when formulated in an ordinary differential equation (ODE) form with linear
time-scaling [27, 28, 29]. In our recent Letter [30], we generalized the GKBA
to quantized bosonic particles. In Ref. [31], we further extended the time-
linear formulation to quantum transport, empowering the Meir-Wingreen formula
[32] for calculating time-dependent currents. The second version of cheers
represents a novel implementation of these ideas adding the time-linear
formulation, the interaction with bosonic particles as well as the coupling to
fermionic baths. The latest improvement of theory [33, 34, 29] was a natural
reason for restructuring of the code and incorporation of new features such as
computation of the molecular integrals based on our earlier implementation
[35]. In what follows we will denote our computational scheme as GKBA+ODE in
order to emphasize the time-linear implementation.

It may be tempting to conclude that GKBA+ODE is a black box method that can be
applied to a variety of systems without any physical insight. One can even
parallel it with another popular approach, time-dependent density functional
theory (TDDFT). Common to them are a small set of observables, linear
propagation time, and reliance on additional approximations. The knowledge of
these approximations distinguishes expert user from a casual practitioner.
They should be carefully selected based on the level of correlations in the
initial and excited state of the studied system [29], and one should be aware
that GKBA, while maintaining the conserving properties of the many-body
approximations, may lead to violations of other physically mandated
properties, such as the positivity of electronic [36] and doublonic [34]
occupation numbers.

The purpose of this contribution is to assist users in taking their first
steps with the cheers code and to explain its structure to developers.
Technical details of the underlying theory and implementation can be found in
the Supporting Information.

##  2 Structure of the code

Mathematically speaking, GKBA+ODE represents an initial value problem (ivp)

$\displaystyle\dot{\bm{\mathit{y}}}(t)$
$\displaystyle=\mathcal{F}\left\\{t,\bm{\mathit{y}}(t)\right\\},$ (1a)
$\displaystyle\bm{\mathit{y}}(t_{0})$ $\displaystyle=\bm{\mathit{y}}_{0},$
(1b)

where $\mathcal{F}$ depends on the time $t$ and is a complicated functional of
the state vector $\bm{\mathit{y}}(t)$. Under the state vector we understand a
collection of electronic, bosonic and embedding correlators that fully
characterize the system under consideration. Its precise composition is not
fixed and depends on the level of approximation.

The goal is to compute the evolution in time of some relevant physical
observables depending on the system state, that is

$\displaystyle\bm{\mathit{O}}(t)$
$\displaystyle=\bm{\mathit{O}}\left\\{\bm{\mathit{y}}(t)\right\\}.$ (2)

Typical examples of observables would be the dipole transition moment,
different contributions to the total energy, electron and doublon occupation
numbers, a set of bosonic observables such as occupation numbers and
displacements, and currents for open systems.

Under some mild conditions, the ivp has a unique solution. Thus, at least
mathematically, the problem is well-posed and requires a mere application of
some well-known methods such as Runge-Kutta (RK) for the approximate solutions
of ordinary differential equations (ODE). While this indeed is the essential
part of the cheers code, the computational complexity of the $\mathcal{F}$
functional grows rapidly with the system size. This is due to the fact that
the state vector $\bm{\mathit{y}}(t)$ contains the one- and two-body density
matrices with dimensions growing as a second and forth power of the system
size, respectively. It represents the major challenge for implementation and
determines the feasibility of the method as well as the maximal propagation
time $t_{\text{f}}$. In what follows, we characterize the system with the
number of electronic $N_{\text{$e$-basis}}$ and bosonic $N_{\text{$b$-basis}}$
basis functions and denote the dimension of the state vector
$\bm{\mathit{y}}(t)$ as $N_{v}$. Already the form of Eqs. (1) and the
dimension $N_{v}$ mandate important design decisions.

###  2.1 General design features

Table 1:  Typical characteristics of organic molecules (amino acids) computed using the pVDZ basis set (for tyrosine 6311G is used). $n_{o}$ and $n_{v}$ represents typical number of occupied and virtual states, respectively. System | Structure | Formula | AO basis | $n_{o}$ | $n_{v}$
---|---|---|---|---|---
Alanine |  |  C3NH7O2 | 125 | 18 | 17
Glycine |  |  C2NH5O2 | 100 | 15 | 15
Tyrosine |  |  C9NH11O3 | 202 | 35 | 35
Proline |  |  C5NH9O2 | 165 | 23 | 24
Phenylalanine |  |  C9NH11O2 | 235 | 32 | 33
Cysteine |  |  C3NH7SO2 | 144 | 21 | 24
Figure 1: Workflow diagram of the GKBA+ODE calculations of the nonequilibrium
dynamics in _molecular systems_ (one can think of them as zero dimensional
systems, and this is reflected in suffixes _0d in the names of respective
modules). Pieces of code belonging to the cheers package are gray-shaded.
Input parameters are framed in blue color, output observables are framed in
red color. Most computationally intensive subroutines are rhs belonging to the
modules Mhf, M2ba, Mgw, Mtph, Mtpp, and Mfad. They all have access to the
module Mgkba_E_R, which contains all relevant system parameters needed for the
time-propagation. Mrk4 is our implementation of the RK4 algorithms. It
contains a stepper rk4step and a driver solveODE, which loops over the times,
estimates propagation error and calls a subroutine for the computation of
observables. Due to its simplicity and efficiency rk4step is not parallelized,
whereas rhs are perfect candidates for parallelization.  Figure 2: Workflow
diagram of the GKBA+ODE calculations of the nonequilibrium dynamics in _model
lattice systems_. Notations as in Fig. 1. Additionally we have a module for
storing bosonic variables Mgkba_B_R and computational bosonic modules
Mehrenfest_0d and Mgd_0d. Notice that electronic and bosonic modules are not
completely independent. For instance Mgd_0d needs to fetch effective
electronic one-body Hamiltonian stored in Mgkba_E_R.

We focus on the scenario where several copies of the state vector
$\bm{\mathit{y}}(t)$ can be stored in the main memory of a computer or be
distributed to each node of a computer cluster. This is by no means obvious:
typical DFT implementations completely avoid the storage of quantities scaling
as a forth power of the system size (such as Coulomb integrals). They operate
with one-body densities with quadratic scaling, whereas the Coulomb matrix
elements in atomic orbital (AO) basis are computed “on the fly”, allowing one
to perform self-consistent calculations even for proteins [37].

This approach is not feasible for a KBE+GKBA theory, which is often formulated
in the molecular orbital (MO) basis making the re-computation of the
Hamiltonian matrix elements at each time-step infeasible. Intricacies of the
basis selection and the basis transformation are discussed by [21], where it
is shown how to deal with a continuum of unbound states that participate in
the photoionization processes. This constitutes a separate preparatory step.
Of course, the fact that
$N_{v}=\mathcal{O}\left(N_{\text{$e$-basis}}^{4}\right)$ puts strong
restrictions of the feasible system sizes. Let us consider a few
representative examples (Tab. 1). As can be seen from the table, the number of
atomic basis functions in the full-electron calculations is significantly
larger than the number of occupied and virtual MO states used in the GKBA+ODE
dynamics. This is achieved by excluding the core and the unoccupied states
with large positive energy. This restriction is not mathematically but
physically motivated considering typical optical excitations using an XUV
laser. For instance, in [29] we performed calculations for the Glycine
molecule using only 5 relevant valence states out of 15 available. This
selection was based on the insight provided by experimental measurements [38]
and theoretical calculations [39, 40].

As a rule, in our calculations we use experimental molecular geometries from
the PubChem [41] online database and Gaussian basis sets from the Basis Set
Exchange [42] online resource. We note, however, that cheers is not restricted
to this scenario, it is conceptualized to work with any programs that can
generate MO integrals (Fig. 1) as well as with model systems (Fig. 2). In
order to get familiar with the second possibility, in Tab. 2 we present
typical resources in order to study the electron localization [43] in the
Hubbard-Holstein model. A novel feature of these calculations is that both
$e$-$e$ and $e$-ph interactions are taken into account.

Table 2:  Largest systems and computational resources (Xeon Gold 5218 CPU @ 2.30GHz) required to study the dynamics of the Holstein-Hubbard model with $N$ sites and physical propagation time $t_{\text{f}}$. We compare timings and sizes of the state vector for a set of electronic (HF — Hartree-Fock, 2B — second Born, $T^{pp}$ — $T$-matrix in the particle-particle channel) and bosonic (Ehrenfest — mean-field, $GD$ — Fan-Migdal approximation) methods.  System | Correlations | State | Time | CPU hours
---|---|---|---|---
$N$ |  $e$-$e$ |  $e$-ph | vector | $t_{\text{f}}$ | $e$-$e$ | $e$-ph
151 | HF | Ehrenfest | 23 254 | 40 | 0.02 | 0.07
151 | HF | $GD$ | 7 046 264 | 40 | 3.8 | 2.0
151 | 2B | $GD$ | 526 954 666 | 40 | 181.2 | 2.0
151 | $T^{pp}$ | Ehrenfest | 519 931 656 | 40 | 482.8 | 0.04
151 | $T^{pp}$ | $GD$ | 526 954 666 | 40 | 488.0 | 2.6

It can be concluded on the basis of these two examples that GKBA+ODE
calculations are routinely possible for amino acid molecules and correlated 1D
clusters. From the structure of the code depicted in Figs. 1, 2, the
separation of the ODE solver from the data generation and the RHS evaluation
($\mathcal{F}\left\\{t,\bm{\mathit{y}}(t)\right\\}$) represents an optimal
architecture amenable to different parallelizations. It is anticipated that
our code can be applied to a variety of physical scenarios, including
ultrafast charge migration in XUV photoexcited halogenated amino acids [44],
2D and 3D model systems. In the next sections we consider some representative
blocks of code according to the workflow diagrams. More detailed description
of the full structure is deferred to the Supporting Information.

###  2.2 Ground state properties and the initialization of the ODE solver

Let us start with molecular systems. As can be seen from the workflow in Fig.
1, a ground state calculation is required in order to generate matrix elements
for subsequent propagation. They are stored in a binary format in the file
with an extension .me containing the following records:

me1p

One-particle matrix elements. They typically comprise the kinetic energy
operator, electron-nuclear repulsion and the mean-field interaction of valence
electrons with the core electrons. Depending on the method, core electrons may
or may not be included in the ground state calculation. The active space for
dynamical calculations is determined by physical considerations.

medp

Dipole transition matrix elements along three Cartesian axes.

me2p

Two-particle matrix elements representing the Coulomb repulsion between the
electrons.

Matrix elements are assumed to be real and they are stored in a _packed form_
in accordance with their symmetries. One of the first tasks is to read the
matrix elements (me1p, medp, me2p) and the system dimensions (the number of
electronic states NP and highest occupied molecular orbital NF) and store them
in the Msys_0d module. We typically assume that HF Hamiltonian is diagonal in
MO basis. This is dictated by practical rather than physical considerations:
for instance one may study correlated molecular dynamics starting from a
weakly correlated ground state in which the initial density matrix is diagonal
and stationary.

However, cheers is not limited to molecular systems. Correspondingly, there
are several Msys_ modules. Information on the _lattice systems_ is stored in
the Msys_lattice module, whereas for restricted Hartree-Fock calculations in
_atomic basis_ , the eri_cheers.x program (part of cheers) makes use of the
Msys_molecule module to store the information about the AO-basis.

The second input file, which is mandatory for all GKBA+ODE calculations, is a
parameter file .par. It consists of several namelists reflecting different
aspects of a simulation:

io

provides names of the matrix elements file and output files for observables.

dynamics

specifies the propagation time-interval, time-step, and the level of treatment
of electronic and bosonic correlations. Currently, only a fixed time-step ODE
solver is used. Although it may not be the fastest option, this choice ensures
that calculations remain reproducible and directly comparable to other codes,
thereby simplifying the debugging process.

sys

gives MO basis size, index of highest occupied molecular orbital, and
dimensionality, number of orbitals per site, number of leads for lattice
systems. The number of read MO states (`np`) may in general be different
(smaller or equal) from the number of MO states in the `.me` file. Also, `np`
may in general be different (greater or equal) from the size of electronic
basis for the GKBA+ODE calculations, $N_{\text{$e$-basis}}$ (`ns`). This gives
one additional flexibility to experiment with different sets of active
molecular orbitals and to exclude irrelevant valence states from the
consideration (known in advance to be always occupied or always empty).

electrons

provides details of electron hopping (nearest and next-nearest) and Coulomb
interaction, periodicity, spinless option, and extension for lattice systems.
The resulting electronic quantities are then stored in Mgkba_E_R. This is one
of the most important modules as it contains all the data for electronic
GKBA+ODE calculations. As is depicted in Fig. 1, the computational routines in
Mhf, M2ba, Mgw, Mtph, Mtpp, and Mfad make use of the matrix elements stored in
Mgkba_E_R and also use the module for intermediate quantities.

phonons

gives phononic frequencies and parameters of $e$-ph interaction. Based on
them, corresponding matrix elements are constructed and stored in the
Mgkba_B_R module. The computational routines in the bosonic modules
Mehrenfest_0d and Mgd_0d make use of the matrix elements and also use the
module for the storage of intermediate quantities.

optical_field

contains parameters of the laser pulses, such as duration, frequency and
shape.

leads

contains tunneling rate matrices for the wide band limit approximation (WBLA),
or details of the leads and tunneling Hamiltonians;

advanced

specifies additional options for the methods used, such as, inclusion of
exchange, freezing of mean field Hamiltonians, presence of acoustic phonons,
treatment of Coulomb integrals with two indices or as complex quantities. Each
of the methods may have different options, all the possibilities are described
in Moptions.

###  2.3 Organization of the state vector

Figure 3: Organization of the state vector. Electronic, bosonic and leads part
with dimensions dim_e, dim_b, and dim_l, respectively, are indicated with
different colors. In the electronic part the electron density matrix
$\rho^{<}$ (Eq. 15) and two-electron correlator $\mathcal{G}^{e}$ (Eq. 17) can
be accessed at off-sets one_body_start and two_body_start, respectively. In
the bosonic part, the bosonic density matrix $\gamma^{<}$ (Eq. 18), electron-
boson correlator $\mathcal{G}^{b}$ (Eq. 21), and expectation values of bosonic
field operators $\phi$ (Eq. 22) can be accessed at off-sets one_body_start,
two_body_start and phi_start, respectively. Finally, the embedding correlator
$\mathcal{G}^{\text{em}}$ (Eq. 23) can be found at the off-set gleads_start.

The code comprises several executables (such as shown in Figs. 1, 2).
Therefore, the state vector $\bm{\mathit{y}}$ is physically located in
different modules Mctrl_e_0d (molecular systems), Mctrl_lattice (lattice
systems) and in Mctrl_transport (lattice systems with embedding). It always
contains the electronic component, and can optionally contain the phononic (or
bosonic) and embedding (or leads) parts, i.e.,

$\displaystyle\bm{\mathit{y}}(t)$
$\displaystyle=\\{\bm{\mathit{y}}_{e}(t),\bm{\mathit{y}}_{b}(t),\bm{\mathit{y}}_{l}(t)\\},$
(3)

as depicted in Fig. 3. In order to properly define each component, let us
consider the most general system Hamiltonian.

#### System Hamiltonian

$\hat{H}(t)=\hat{H}_{\text{el}}(t)+\hat{H}_{\text{bos}}(t)+\hat{H}_{\text{el-
bos}}(t)+\hat{H}_{\text{leads}}+\hat{H}_{\text{tunnel}}.$ (4)

It contains the electronic part

$\displaystyle\hat{H}_{\text{el}}(t)$
$\displaystyle=\sum_{ij}h_{ij}(t)\hat{d}_{i}^{\dagger}\hat{d}_{j}+\frac{1}{2}\sum_{ijmn}v_{ijmn}(t)\hat{d}_{i}^{\dagger}\hat{d}_{j}^{\dagger}\hat{d}_{m}\hat{d}_{n},$
(5)

comprising a one-body term ($h^{\dagger}=h$) accounting for the kinetic energy
as well as the interaction with nuclei and possible external fields. The time-
dependence of $h_{ij}(t)=h^{0}_{ij}+\langle i|V_{\rm
gate}(t)+\hat{\mathbf{d}}\cdot{\mathbf{E}_{\rm pulse}}(t)|j\rangle$ is due to
a time-dependent gate voltage $V_{\mathrm{gate}}(t)$ and to a possible laser
pulse $\mathbf{E}_{\rm pulse}(t)$ coupled to the electronic dipole operator
$\hat{\mathbf{d}}$. Furthermore, there is a two-body term accounting for the
Coulomb interaction between the electrons. The time-dependence of the Coulomb
matrix elements $v_{ijmn}(t)=v_{ijmn}s(t)$ could be due to the adiabatic
switching protocol $s(t)$ adopted to generate a correlated initial state.
Latin indices denote one-electron states comprising orbital and spin degrees
of freedom, i. e., $i=(\bm{i},\sigma)$. As can be seen from Figs. 1, 2, $h$
and $v$ matrix elements are stored in the Mgkba_E_R module.

We write the bosonic Hamiltonian as

$\displaystyle\hat{H}_{\text{bos}}(t)=\sum_{\mu\nu}\Omega_{\mu\nu}(t)\hat{\phi}_{\mu}\hat{\phi}_{\nu},$
(6)

where $\Omega^{\dagger}=\Omega$ may depend on time, such as in the parametric
phonon drivings [45] setup. The annihilation and creation operators for a
bosonic mode $\mu$, i.e.,
$\hat{a}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}$
and
$\hat{a}^{{\dagger}}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}$,
are arranged into a vector
$(\hat{x}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}},\hat{p}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}})$
where
$\hat{x}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}=(\hat{a}^{{\dagger}}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}+\hat{a}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}})/\sqrt{2}$
are the position operators and
$\hat{p}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}=i(\hat{a}^{{\dagger}}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}-\hat{a}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}})/\sqrt{2}$
are the momentum operators. The greek index $\mu=(\mbox{\boldmath$\mu$},\xi)$
is then used to specify the bosonic mode and the component of the vector:
$\hat{\phi}_{\mu}=\hat{x}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}$
for $\xi=1$ and
$\hat{\phi}_{\mu}=\hat{p}_{\mbox{\scalebox{0.7}{$\scriptscriptstyle\mbox{\boldmath$\mu$}$}}}$
for $\xi=2$. The electronic and bosonic subsystems interact through

$\displaystyle\hat{H}_{\text{el-bos}}(t)$
$\displaystyle=\sum_{\mu,ij}g_{\mu,ij}(t)\hat{d}_{i}^{\dagger}\hat{d}_{j}\hat{\phi}_{\mu};$
(7)

therefore electrons can be coupled to both the mode coordinates and momenta.
We allow $g$ to depend on time for possible adiabatic switchings, i. e.,
$g_{\mu,ij}(t)=g_{\mu,ij}s(t)$. $\Omega$ and $g$ matrix elements are stored in
the Mgkba_B_R module.

Similarly to the bosonic part, the lead part of the total Hamiltonian is
optional: not every calculation is performed with transport setup in mind,
moreover the leads and the tunneling Hamiltonian parts do not explicitly enter
the formalism in the wide-band limit approximation. Nonetheless, we provide
some possible explicit forms of these parts in order to give a microscopic
definition of the tunneling-rate matrices. Thus, the leads Hamiltonian can be
parametrized as

$\displaystyle\hat{H}_{\text{leads}}$
$\displaystyle=\sum_{k,\alpha,\sigma}\epsilon_{k\alpha}\hat{c}_{k\alpha,\sigma}^{\dagger}\hat{c}_{k\alpha,\sigma}.$
(8)

At the moment we assume one-dimensional momenta $k$. Here $\alpha$ denotes the
lead number and $\sigma$ is the spin projection. For a one-dimensional tight-
binding chain (with on-site energy $a_{\alpha}$ and hopping $b_{\alpha}$) the
electron energy dispersion $\epsilon_{k\alpha}$ is parametrized as

$\displaystyle\epsilon_{k\alpha,\sigma}$
$\displaystyle=a_{\alpha}+2|b_{\alpha}|\cos\left(\frac{k\pi}{N_{k}+1}\right),$
$\displaystyle 1$ $\displaystyle\leq k\leq N_{k}.$ (9)

As can be seen, $a$ and $b$ are lead-dependent, whereas we assume the same
number of momenta $N_{k}$ per each lead. Here, leads are not spin-polarized,
but it would be straightforward to consider, e. g., ferromagnetic leads via
modified lead energies. Upon the application of a time-dependent bias voltage
$V_{\alpha}(t)$, the lead energies are modified as

$\displaystyle\tilde{\epsilon}_{k\alpha,\sigma}(t)$
$\displaystyle=\epsilon_{k\alpha,\sigma}+V_{\alpha}(t).$ (10)

The tunneling Hamiltonian can be parametrized as:

$\displaystyle\hat{H}_{\text{tunnel}}(t)$
$\displaystyle=\sum_{k,\alpha}\sum_{\bm{i},\sigma}T_{\bm{i},k\alpha}(t)\left[\hat{d}_{\bm{i},\sigma}^{\dagger}\hat{c}_{k\alpha,\sigma}+h.c.\right].$
(11)

The coupling strengths are in general time- and lead-dependent, but it is
assumed that for each lead $\alpha$ there is a _single contact site_ , i.e.,

$\displaystyle T_{\bm{i},k\alpha}(t)$
$\displaystyle=t_{\alpha}\delta_{\bm{i}C_{\alpha}}\sqrt{\frac{2}{N_{k}+1}}\sin\left(\frac{k\pi}{N_{k}+1}\right)s_{\alpha}(t).$
(12)

$C_{\alpha}$ represents the quantum number (such as site number and band
index) of the electronic state that couples to the lead $\alpha$.
$s_{\alpha}(t)$ is the switch-on function between the system and the lead
$\alpha$, see Fig. 2. Next we define the tunneling-rate matrix

$\displaystyle\Gamma_{\alpha,\bm{ij}}(\omega;t)$
$\displaystyle=2\pi\sum_{k=1}^{N_{k}}T_{\bm{i},k\alpha}(t)\delta(\omega-\epsilon_{k\alpha})T_{\bm{j},k\alpha}^{*}(t),$
(13)

which is in general frequency-dependent, but in WBLA it reduces to

$\displaystyle\Gamma_{\alpha,\bm{ij}}(\omega;t)$
$\displaystyle=\gamma_{\alpha}\delta_{\bm{i}C_{\alpha}}\delta_{\bm{j}C_{\alpha}}s_{\alpha}^{2}(t).$
(14)

#### Parts of the state vector

The electronic component $\bm{\mathit{y}}_{e}$ contains i) the electron
density matrix

$\displaystyle\rho_{ij}^{<}$ $\displaystyle=-iG^{<}_{ij}(t,t),$ (15)
$\displaystyle G^{<}_{ij}(t,t^{\prime})$
$\displaystyle=i\langle\hat{d}^{\dagger}_{j}(t^{\prime})\hat{d}_{i}(t)\rangle,$
(16)

proportional to the equal-time lesser Green’s function $G^{<}$; and ii) a two-
electron correlator

$\displaystyle\mathcal{G}^{e}_{imjn}(t)$
$\displaystyle=-\langle\hat{d}^{{\dagger}}_{n}(t)\hat{d}^{{\dagger}}_{j}(t)\hat{d}_{i}(t)\hat{d}_{m}(t)\rangle_{c},$
(17)

The subscript “$c$” in the averages signifies that only the correlated part
must be retained.

The bosonic component $\bm{\mathit{y}}_{b}$ contains i) the bosonic density
matrix $\gamma^{<}_{\mu\nu}$

$\displaystyle\gamma_{\mu\nu}^{<}(t)$ $\displaystyle=iD^{<}_{\mu\nu}(t,t),$
(18) $\displaystyle D^{<}_{\mu\nu}(t,t^{\prime})$
$\displaystyle=D^{>}_{\nu\mu}(t^{\prime},t)=-i\langle\Delta\hat{\phi}_{\nu}(t^{\prime})\Delta\hat{\phi}_{\mu}(t)\rangle,$
(19) $\displaystyle\Delta\hat{\phi}_{\nu}(t)$
$\displaystyle\equiv\hat{\phi}_{\nu}(t)-\langle\hat{\phi}_{\nu}(t)\rangle,$
(20)

proportional to the equal-time lesser bosonic Green’s function
$D^{<}_{\mu\nu}$; ii) a higher-order electron-boson correlator

$\displaystyle\mathcal{G}^{b}_{\mu,ij}(t)$
$\displaystyle=\langle\hat{d}^{{\dagger}}_{j}(t)\hat{d}_{i}(t)\hat{\phi}_{\mu}(t)\rangle_{c},$
(21)

and iii) the expectation values of the bosonic field operators:

$\displaystyle\phi_{\nu}(t)$
$\displaystyle\equiv\langle\hat{\phi}_{\nu}(t)\rangle.$ (22)

Matrices $\gamma^{<}$, $D^{<}$ have dimensions $N_{\text{$b$-basis}}\times
N_{\text{$b$-basis}}$.

To deal with open systems, the leads part $\bm{\mathit{y}}_{l}$ contains the
embedding correlator $\mathcal{G}^{\text{em}}$, which can be completely
specified in terms of the leads and tunneling Hamiltonians, together with the
quasi-particle propagators of the coupled system. This is due to the leads
being non-interacting and in contrast to other higher-order correlators such
as Eq. (17). Here, we consider the wide-band limit approximation by expressing
the embedding correlator within the so-called pole expansion scheme [31],

$\displaystyle\mathcal{G}_{l\alpha}^{\text{em}}(t)=\int
d\bar{t}s_{\alpha}(\bar{t}){\rm e}^{-i\phi_{\alpha}(t,\bar{t})}{\rm
e}^{-i\mu(t-\bar{t})}{\rm
e}^{-\zeta_{l}(t-\bar{t})/\beta}G^{\text{A}}(\bar{t},t).$ (23)

While we will next specify all the components of Eq. (23), we refer the Reader
to the Supporting Information and our recent letter [31] for additional
details. The embedding correlator, besides depending on the lead index
$1\leq\alpha\leq N_{\text{leads}}$ contains the dependence on the pole index
$1\leq l\leq N_{p}$ in the expansion of the Fermi distribution function

$\displaystyle f(x)=\frac{1}{{\rm e}^{\beta x}+1}$
$\displaystyle=\frac{1}{2}-\lim_{N_{p}\to\infty}\sum_{l=1}^{N_{p}}\eta_{l}\left(\frac{1}{\beta
x+i\zeta_{l}}+\frac{1}{\beta x-i\zeta_{l}}\right),$ (24)

where $\eta$ and $\pm i\zeta$ are the residues and poles ($\zeta>0$),
respectively. We typically use the expansion coefficients generated by solving
an eigenvalue problem of a specific, tridiagonal matrix [46]. Additionally
$\phi_{\alpha}(t,t^{\prime})$ is the bias-voltage phase factor

$\displaystyle\phi_{\alpha}(t,t^{\prime})$
$\displaystyle\equiv\int_{t^{\prime}}^{t}d\bar{t}V_{\alpha}(\bar{t}),$ (25)

$\mu$ is the chemical potential, and

$\displaystyle G^{\text{R}/\text{A}}(t,t^{\prime})=\mp
i\theta[\pm(t-t^{\prime})]\mathrm{T}{\rm
e}^{-i\int_{t^{\prime}}^{t}d\bar{t}h_{\text{eff}}^{e}(\bar{t})}$ (26)

are the quasiparticle retarded/advanced propagators of the coupled system.
Notice, that at variance with closed systems, the effective Hamiltonian is in
general non-self-adjoint

$\displaystyle h_{\text{eff}}^{e}(t)$
$\displaystyle=h^{e}(t)-\frac{i}{2}\sum_{\alpha}\Gamma_{\alpha}s_{\alpha}^{2}(t)=h^{e}(t)-\frac{i}{2}\Gamma(t).$
(27)

Matrices $h^{e}$, $\Gamma$, $G^{\text{R}/\text{A}}$,
$\mathcal{G}_{l\alpha}^{\text{em}}$ have dimensions
$N_{\text{$e$-basis}}\times N_{\text{$e$-basis}}$. We conclude this section by
defining the mean-field electronic Hamiltonian

$\displaystyle h^{e}_{ij}(t)$ $\displaystyle=h_{ij}(t)+V^{\rm HF}_{ij}(t),$
(28) $\displaystyle V^{\rm HF}_{ij}(t)$
$\displaystyle=\sum_{mn}w_{imnj}\rho^{<}_{nm}(t),$ (29)

with $w_{imnj}=v_{imnj}-v_{imjn}$ being the antisymmetrized interaction. The
storage space for $\bm{\mathit{y}}$ — $\mathcal{O}(N_{v})$ — is dominated by
the $\mathcal{G}^{e}$ correlator — $\mathcal{O}(N_{\text{$e$-basis}}^{4})$,
the electron-boson correlator $\mathcal{G}^{b}$ —
$\mathcal{O}(N_{\text{$e$-basis}}^{2}N_{\text{$b$-basis}})$, and the embedding
correlator — $\mathcal{O}(N_{\text{$e$-basis}}^{2}N_{\text{leads}}N_{p})$,
where $N_{\text{$b$-basis}}$ is double of the number of bosonic modes, and the
number of poles in the expansion of the Fermi distribution function is
typically below 50.

###  2.4 Computational subroutines

In the previous section we introduced the electronic, bosonic and leads parts
of the state vector. Equations of motion covering the combinations of
electronic-bosonic and electronic-leads degrees of freedom were derived in
Ref. [33] and Ref. [31], respectively. These equations demonstrate the
conserving properties [30] and symmetries of underlying theories. It follows,
for instance, that density matrices $\rho^{<}$ (Eq. 15) and $\gamma^{<}$ (Eq.
18) are self-adjoint and that higher-order correlators can be likewise written
in matrix form, by combing electronic indices in pairs. From the electronic
correlator tensor $\mathcal{G}^{e}_{imjn}$ (Eq. 17) a matrix
$\bm{\mathcal{G}}^{e}(t)$ can be constructed in three different ways depending
on the approximation used: $GW$, $T^{ph}$, or $T^{pp}$. These approximations
contain an infinite sequence of diagrams corresponding to the solution of the
Bethe-Salpeter equation in three different channels. $GW$ approximation
corresponds to the treatment of $\overline{ph}$ channel, $T^{pp}$
approximation — $pp$ channel, and $T^{ph}$ — $ph$ channel, where we refer to
the naming convention used in the parquet theory [47]. Additionally, it is
possible to take even more complicated diagrams into account by including
exchange contributions [33]. In every such form the matrix
$\bm{\mathcal{G}}^{e}(t)$ can be shown to be self-adjoint [29]. This matrix
form of the GKBA+ODE method leads to important simplifications of the
resulting equations, for instance only the upper triangular parts must be
computed. For the explicit form of the equations we refer to corresponding
papers [33, 31] and focus here more on their general functional form.

The evaluation of $\mathcal{F}\left\\{t,\bm{\mathit{y}}(t)\right\\}$ in Eq.
(1a), which constitutes the main computational task, can be partitioned as

$\displaystyle\mathcal{F}\left\\{t,\bm{\mathit{y}}(t)\right\\}$
$\displaystyle=\left\\{\mathcal{F}_{e}\left\\{t,\bm{\mathit{y}}_{e}(t),\bm{\mathit{y}}_{b}(t),\bm{\mathit{y}}_{l}(t)\right\\},\mathcal{F}_{b}\left\\{t,\bm{\mathit{y}}_{e}(t),\bm{\mathit{y}}_{b}(t)\right\\},\mathcal{F}_{l}\left\\{t,\bm{\mathit{y}}_{e}(t),\bm{\mathit{y}}_{l}(t)\right\\}\right\\}.$
(30)

As can be seen from this form, the evaluation of one part requires the
knowledge of other parts. This structure is present even at the mean-field
level, where the equation of motion for the bosonic and leads correlators
depends on the electron density matrix. In correlated cases, there are
feedbacks of these degrees of freedom on electrons, and the electron density
matrix $\rho^{<}$ enters the evaluation of the higher-order correlators. These
complicated relations are shown as arrows in Figs. 1, 2. The subroutines
implementing different electronic, bosonic and embedding approximations are
grouped in modules such as:

Mhf_0d

mean-field for electrons,

M2ba_0d

second Born approximation for electrons,

Mgw_0d and Mgw_hub

$GW$ approximation for molecular and lattice systems,

Mtpp_0d and Mtpp_hub

$T^{pp}$ approximation for molecular and lattice systems,

Mtph_0d and Mtph_hub

$T^{ph}$ approximation for molecular and lattice systems,

Mfad_0d

three-particle Faddeev method for molecular systems,

Mehrenfest_0d and Mgd_0d

Ehrenfest and Fan-Migdal (or $GD$) approximation for phonons,

Membedding

the spectral decomposition and the pole expansion methods for leads.

In addition to the state vector, the rhs subroutines (computing right hand
sides of the equations of motion) in these modules rely on the common data
(matrix elements) stored in Mgkba_E_R and Mgkba_B_R modules. The code is
organized in such a way that subroutines implementing different approximations
are completely interchangeable. For instance, every evaluation starts with a
transformation between the state vector (real) $\bm{\mathit{y}}$ and the
matrix form (complex) of $\rho^{<}$, $\bm{\mathcal{G}}^{e}(t)$ for electrons,
$\gamma^{<}$, $\bm{\mathcal{G}}^{b}(t)$, $\phi$ for bosons, and
$\mathcal{G}^{\text{em}}$ for leads. This computational overhead is needed for
three reasons: i) simplicity and universality of the differential equation
solver solveODE; ii) speed-up due to automatic loops unrolling in the Runge-
Kutta algorithm, rk4step; iii) the computational complexity of the
transformation operations is well below time-critical operations, which scale
as fifth or higher power of the system size. There is also small memory saving
due to the use of symmetry, which may play an important role in distributed
computing. Therefore, it is possible to achieve a great variety of methods
without code duplication. The details of this design are presented in
Supporting Information. Let us now demonstrate the interplay of electronic
correlations, electron-phonon interaction and tunneling to the leads in a
model calculation.

##  3 Illustrative example

We consider here a 1D Hubbard model at low filling interacting at every site
with a phonon mode. This is also known as the Holstein-Hubbard model, which is
often studied in the context of polaron formation and electron localization
[43] and field-driven heating [48]. We add here another aspect to these
studies, namely the presence of leads connected to the two terminal sites
($N_{\text{leads}}=2$). This designates our model as a prototype for a
photovoltaic device encompassing the phenomena of creation of nonequilibrium
carriers, formation of the polaronic quasiparticles, and charge transport
through the active material to leads. We focus here only on modelling of the
conduction band. A realistic model of the whole photovoltaic device should
also comprise the valence band, from where the electrons are excited, and the
optical pulse. This calculation is fully within the capability of the code.
However, here we deal with a simpler one-band scenario and perform
calculations for shorter chains as compared to previous studies [34] — only 35
sites ($N_{v}=1792350$) — but they capture all relevant physical phenomena. As
a benefit, all the calculations can be performed on a typical laptop within
five minutes time. The Hamiltonian is given by

$\displaystyle\hat{H}_{\text{h-h}}$
$\displaystyle=-t_{\text{nn}}\sum_{\sigma=\uparrow,\downarrow}\sum_{<\bm{i},\bm{j}>}\hat{d}_{\bm{i}\sigma}^{\dagger}\hat{d}_{\bm{j}\sigma}+U\sum_{\bm{i}}\hat{n}_{\bm{i}\uparrow}\hat{n}_{\bm{i}\downarrow}+\sum_{\bm{i}}\left\\{\Omega\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{i}}+g(\hat{a}^{\dagger}_{\bm{i}}+\hat{a}_{\bm{i}})\hat{n}_{\bm{i}}\right\\},$
(31)

where $\hat{n}_{\bm{i}}=\sum_{\sigma}\hat{n}_{\bm{i}\sigma}$,
$\hat{n}_{\bm{i}\sigma}=\hat{d}_{\bm{i}\sigma}^{\dagger}\hat{d}_{\bm{i}\sigma}$,
and we set $t_{\text{nn}}=\Omega=1$, and $U=2$, $g=0.5$. The dynamics is
triggered by the creation of an electron pair at the central lattice site.
Initially the lattice is empty. In the absence of leads the electronic wave
packet spreads on the lattice, whereby the electron group velocity is at most
$2t_{\text{nn}}$. Two electrons occupying the same lattice site form a
doublon, a long-living excitation, the decay of which is suppressed by energy
conservation. Its group velocity is in general different from the electronic
one and in the strongly interacting limit can be estimated as
$4t_{\text{nn}}^{2}/U$. Propagation of electrons and doublons at low filling
can be accurately described by the $T^{pp}$ approximation [34]. Both
velocities are renormalized by the $e$-ph interaction.

If such a system is subject to open boundary conditions, the electronic wave-
packet will hit the chain boundaries and reflect after propagating for a
certain time. If the system is connected to leads, the charge will tunnel into
or from the central region, and the number of electrons in the system will not
be preserved. In the present scenario, we study how $e$-$e$ and $e$-ph
interactions influence the charge flow to leads in the wide band limit
approximation. The chemical potential of the leads and central system is
assumed to be equal, however, we set the lead voltage to a large negative
value ($V_{\alpha}=-5$). This enables electrons to leave the chain. The
dynamics of the system can be understood based on the analysis of total
populations in Fig. 4.

Figure 4: Total number of electrons and doublons in a 35-sites system
connected to leads with $V_{\alpha}=-5$, $\gamma_{\alpha}=1$, and $\beta=100$.
Different line-styles denote different levels of treatment of $e$-$e$ and
$e$-ph interactions. Electrons are characterized by the nearest-neighbor
hopping $t_{\text{nn}}=1$ and on-site Hubbard interaction $U=2$. The phonon
frequency is set to $\Omega=1$, and the electron-phonon coupling constant for
calculations within the Fan-Migdal approach (blue solid lines) or the
Ehrenfest approach (red dash-dotted line) is set to $g=0.5$. The orange dashed
line denotes results without the electron-phonon interaction ($g=0$).
Electrons are treated either at the Hartree-Fock (HF) or the $T$-matrix in the
particle-particle channel level of approximation.

The considered system is not spin polarized. Therefore, only the spatial
components of $\rho^{<}$ and $\mathcal{G}^{e}$ are propagated. They are
defined as

$\displaystyle\rho_{\bm{\mathit{i}}\bm{\mathit{j}}}^{<}$
$\displaystyle\equiv\rho_{\bm{\mathit{i}}\sigma,\bm{\mathit{j}}\sigma}^{<},$
(32)
$\displaystyle\mathcal{G}^{e}_{\bm{\mathit{i}}\bm{\mathit{m}}\bm{\mathit{j}}\bm{\mathit{n}}}$
$\displaystyle\equiv\mathcal{G}^{e}_{\bm{\mathit{i}}\sigma,\bm{\mathit{m}}\sigma,\bm{\mathit{j}}\sigma,\bm{\mathit{n}}\sigma},$
(33)

where $\sigma=\uparrow,\downarrow$. The electronic and doublonic populations
are defined then as respective diagonal elements

$\displaystyle n_{\bm{\mathit{i}}}$
$\displaystyle\equiv\langle\hat{n}_{\bm{\mathit{i}},\sigma}\rangle=\rho_{\bm{\mathit{i}}\bm{\mathit{i}}}^{<},$
(34) $\displaystyle d_{\bm{\mathit{i}}}$
$\displaystyle\equiv\langle\hat{n}_{\bm{\mathit{i}},\uparrow}\hat{n}_{\bm{\mathit{i}},\downarrow}\rangle-
n_{\bm{\mathit{i}}}^{2}=-\mathcal{G}^{e}_{\bm{\mathit{i}}\bm{\mathit{i}}\bm{\mathit{i}}\bm{\mathit{i}}}.$
(35)

Initially the total number of spin-up or spin-down electrons is

$\displaystyle n_{e}$
$\displaystyle=\sum_{\bm{\mathit{i}}}n_{\bm{\mathit{i}}}=1,$ (36)

and the system is uncorrelated, therefore the number of doublons is zero

$\displaystyle n_{d}$
$\displaystyle=\sum_{\bm{\mathit{i}}}d_{\bm{\mathit{i}}}=0,$ (37)

notice that uncorrelated part is subtracted in the definition of
$d_{\bm{\mathit{i}}}$.

It is assumed that leads with $\gamma_{\alpha}=1$ and inverse temperature
$\beta=100$ are attached suddenly at $t=0$, i. e., $s_{\alpha}(t)=\theta(t)$,
as the propagation starts at $t_{i}=0$. Since the central system and the leads
are initially not equilibrated, the relaxation gives rise to a small flow of
charge in and out of the system at $t\simeq 0.5$. After this kink, $n_{e}$
remains constant, and the build up of $e$-$e$ and $e$-ph correlations takes
place. This is reflected in the rapid growth of the doublon population,
primarily driven by pure electronic scattering, and its oscillations resulting
from interactions with phonons. Simultaneously, electrons and doublons
propagate through the lattice with distinct group velocities until they reach
the boundaries. It is important to note that the travel times for electrons
and doublons are significantly different (compare the length of plateau in the
plots of $n_{e}$ and $n_{d}$). After reaching the boundaries, electrons tunnel
to the leads indirectly reducing the number of doublons in the system. The
rate of tunneling strongly depends on the presence of phonons and on the level
of approximation used. It can be seen, for instance, that $T^{pp}$
approximation for the electron and Fan-Migdal (or $GD$) approximation for the
phonon self-energies, respectively, predicts the slowest decline of
populations (smallest current). In the absence of $e$-ph coupling ($g=0$)
electrons leak to the leads faster, whereas the mean-field treatment of both
interactions, which does not capture the formation of bound doublon or polaron
states, results in the fastest depletion of the central region.

##  4 Conclusions

We presented a NEGF toolbox for the simulation of the coupled electron-boson
dynamics in correlated materials including the treatment of open systems with
fermionic bath. The merits of our implementation in the cheers package are as
follows: (1) all fundamental conservation laws are satisfied independently of
the method [30]; (2) the ODE nature of the EOM allows one to address phenomena
occurring at different time scales through a save-and-restart procedure; and
(3) as a byproduct of the calculations we have access to the spatially
nonlocal correlators $\bm{\mathcal{G}}^{e}(t)$ and $\bm{\mathcal{G}}^{b}(t)$,
as well as to the embedding correlator $\mathcal{G}^{\text{em}}$ that allows
evaluation of currents according to the Meir-Wingreen formula [31]. The ODE
formulation naturally accommodates parallel computations, adaptive time-
stepping implementations, and restart protocols, thereby paving the way for
conducting first principles simulations of multiscale systems. As an
illustration of our method, we investigated a 1D Holstein-Hubbard model in a
setup corresponding to the conduction band of a photovoltaic device. We found
that the formation of bound doublonic and polaronic states reduces the photo-
current.

## Supporting Information

Supporting Information is available from the Wiley Online Library or from the
authors. At present, the code is given only for collaborative purposes. The
interested reader can contact the authors of the code.

## Acknowledgements

This research is part of the project No. 2021/43/P/ST3/03293 co-funded by the
National Science Centre and the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement no.
945339 (Y.P.). G.S. and E.P. acknowledge funding from MIUR PRIN Grant No.
20173B72NB, from the INFN17-Nemesys project. G.S. and E.P. acknowledge Tor
Vergata University for financial support through projects ULEXIEX and TESLA.

## References

  * [1] D. N. Basov, R. D. Averitt, D. Hsieh,  _Nat. Mater._ 2017, _16_ , 11 1077.
  * [2] G. Stefanucci, R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ ,  Cambridge University Press, Cambridge, 2013.
  * [3] N.-H. Kwong, M. Bonitz,  _Phys. Rev. Lett._ 2000, _84_ , 8 1768.
  * [4] N. E. Dahlen, R. van Leeuwen,  _Phys. Rev. Lett._ 2007, _98_ , 15 153004.
  * [5] P. Myöhänen, A. Stan, G. Stefanucci, R. van Leeuwen,  _Eurphys. Lett._ 2008, _84_ 67001\.
  * [6] M. Galperin, S. Tretiak,  _J. Chem. Phys._ 2008, _128_ , 12 124705.
  * [7] P. Myöhänen, A. Stan, G. Stefanucci, R. van Leeuwen,  _Phys. Rev. B_ 2009, _80_ , 11 115107.
  * [8] M. P. von Friesen, C. Verdozzi, C.-O. Almbladh,  _Phys. Rev. Lett._ 2009, _103_ , 17 176404.
  * [9] N. Bittner, D. Golež, H. U. R. Strand, M. Eckstein, P. Werner,  _Phys. Rev. B_ 2018, _97_ , 23 235125.
  * [10] M. Schüler, J. Berakdar, Y. Pavlyukh,  _Phys. Rev. B_ 2016, _93_ , 5 054303.
  * [11] M. Schüler, D. Golež, Y. Murakami, N. Bittner, A. Herrmann, H. U. Strand, P. Werner, M. Eckstein,  _Comp. Phys. Commun._ 2020, _257_ 107484\.
  * [12] J. Kaye, D. Golez,  _SciPost Physics_ 2021, _10_ , 4 091.
  * [13] F. Meirinhos, M. Kajan, J. Kroha, T. Bode,  _SciPost Physics Core_ 2022, _5_ , 2 030.
  * [14] X. Dong, E. Gull, H. U. R. Strand,  _Phys. Rev. B_ 2022, _106_ , 12 125153.
  * [15] P. Lipavský, V. Špička, B. Velický,  _Phys. Rev. B_ 1986, _34_ , 10 6933.
  * [16] E. V. n. Boström, A. Mikkelsen, C. Verdozzi, E. Perfetto, G. Stefanucci,  _Nano Lett._ 2018, _18_ , 2 785.
  * [17] F. Covito, E. Perfetto, A. Rubio, G. Stefanucci,  _Phys. Rev. A_ 2018, _97_ , 6 061401(R).
  * [18] S. Latini, E. Perfetto, A.-M. Uimonen, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. B_ 2014, _89_ , 7 075306.
  * [19] E. P. Månsson, S. Latini, F. Covito, V. Wanie, M. Galli, E. Perfetto, G. Stefanucci, H. Hübener, U. De Giovannini, M. C. Castrovilli, et al.,  _Commun. Chem._ 2021, _4_ , 1 73.
  * [20] E. Perfetto, D. Sangalli, A. Marini, G. Stefanucci,  _J. Phys. Chem. Lett._ 2018, _9_ , 6 1353.
  * [21] E. Perfetto, D. Sangalli, M. Palummo, A. Marini, G. Stefanucci,  _J. Chem. Theory Comput._ 2019, _15_ , 8 4526.
  * [22] E. Perfetto, A. Trabattoni, F. Calegari, M. Nisoli, A. Marini, G. Stefanucci,  _J. Phys. Chem. Lett._ 2020, 9.
  * [23] E. Perfetto, A.-M. Uimonen, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. A_ 2015, _92_ , 3 033419.
  * [24] E. Perfetto, Y. Pavlyukh, G. Stefanucci,  _Phys. Rev. Lett._ 2022, _128_ , 1 016801.
  * [25] E. Perfetto, G. Stefanucci,  _Nano Lett._ 2023, _23_ , 15 7029.
  * [26] E. Perfetto, G. Stefanucci,  _J. Phys. Condens. Matter_ 2018, _30_ , 46 465901.
  * [27] N. Schlünzen, J.-P. Joost, M. Bonitz,  _Phys. Rev. Lett._ 2020, _124_ , 7 076601.
  * [28] J.-P. Joost, N. Schlünzen, M. Bonitz,  _Phys. Rev. B_ 2020, _101_ , 24 245101.
  * [29] Y. Pavlyukh, E. Perfetto, G. Stefanucci,  _Phys. Rev. B_ 2021, _104_ , 3 035124.
  * [30] D. Karlsson, R. van Leeuwen, Y. Pavlyukh, E. Perfetto, G. Stefanucci,  _Phys. Rev. Lett._ 2021, _127_ , 3 036402.
  * [31] R. Tuovinen, Y. Pavlyukh, E. Perfetto, G. Stefanucci,  _Phys. Rev. Lett._ 2023, _130_ 246301\.
  * [32] Y. Meir, N. S. Wingreen,  _Phys. Rev. Lett._ 1992, _68_ , 16 2512.
  * [33] Y. Pavlyukh, E. Perfetto, D. Karlsson, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. B_ 2022, _105_ , 12 125134.
  * [34] Y. Pavlyukh, E. Perfetto, D. Karlsson, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. B_ 2022, _105_ , 12 125135.
  * [35] Y. Pavlyukh, J. Berakdar,  _Comp. Phys. Commun._ 2013, _184_ 387\.
  * [36] J.-P. Joost, N. Schlünzen, H. Ohldag, M. Bonitz, F. Lackner, I. Březinová,  _Phys. Rev. B_ 2022, _105_ , 16 165155.
  * [37] D. J. Cole, N. D. M. Hine,  _J. Phys. Condens. Matter_ 2016, _28_ , 39 393001.
  * [38] R. Weinkauf, P. Schanen, A. Metsala, E. Schlag, M. Bürgle, H. Kessler,  _The Journal of Physical Chemistry_ 1996, _100_ , 47 18567.
  * [39] A. I. Kuleff, J. Breidbach, L. S. Cederbaum,  _J. Chem. Phys._ 2005, _123_ , 4 044111.
  * [40] D. Ayuso, A. Palacios, P. Decleva, F. Martín,  _Phys. Chem. Chem. Phys._ 2017, _19_ , 30 19767.
  * [41] S. Kim, J. Chen, T. Cheng, A. Gindulyte, J. He, S. He, Q. Li, B. A. Shoemaker, P. A. Thiessen, B. Yu, L. Zaslavsky, J. Zhang, E. E. Bolton,  _Nucleic Acids Res._ 2021, _49_ , D1 D1388.
  * [42] B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson, T. L. Windus,  _J. Chem. Inf. Model._ 2019, _59_ , 11 4814.
  * [43] B. Kloss, D. R. Reichman, R. Tempelaar,  _Phys. Rev. Lett._ 2019, _123_ , 12 126601.
  * [44] G. Bergamaschi, L. Lascialfari, A. Pizzi, M. I. Martinez Espinoza, N. Demitri, A. Milani, A. Gori, P. Metrangolo,  _Chem. Comm._ 2018, _54_ , 76 10718.
  * [45] Y. Murakami, N. Tsuji, M. Eckstein, P. Werner,  _Phys. Rev. B_ 2017, _96_ 045125\.
  * [46] J. Hu, R.-X. Xu, Y. Yan,  _J. Chem. Phys._ 2010, _133_ 101106\.
  * [47] K. Held,  In E. Pavarini, E. Koch, D. Vollhardt, A. I. Lichtenstein, editors, _DMFT at 25: infinite dimensions: lecture notes of the Autumn School on Correlated Electrons 2014_. Forschungszentrum Jülich, Zentralbibliothek, Verl, Jülich, 2014.
  * [48] M. Weber, J. K. Freericks,  _Phys. Rev. Lett._ 2023, _130_ , 26 266401.