# Performance of wave function and Green’s function methods for non-
equilibrium many-body dynamics

Cian C. Reeves  Department of Physics, University of California, Santa
Barbara, Santa Barbara, CA 93117  Gaurav Harsha  Department of Chemistry,
University of Michigan, Ann Arbor, Michigan 48109, USA  Avijit Shee
Department of Chemistry, University of California, Berkeley, USA  Yuanran Zhu
Applied Mathematics and Computational Research Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Chao Yang  Applied Mathematics
and Computational Research Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  K Birgitta Whaley  Department of Chemistry,
University of California, Berkeley, USA  Berkeley Center for Quantum
Information and Computation, Berkeley  Dominika Zgid  Department of Chemistry,
University of Michigan, Ann Arbor, Michigan 48109, USA  Department of Physics,
University of Michigan, Ann Arbor, Michigan 48109, USA  Vojtěch Vlček
Department of Chemistry and Biochemistry, University of California, Santa
Barbara, Santa Barbara, CA 93117  Department of Materials, University of
California, Santa Barbara, Santa Barbara, CA 93117

###### Abstract

Theoretical descriptions of non-equilibrium dynamics of quantum many-body
systems essentially employ either (i) explicit treatments, relying on
truncation of the expansion of the many-body wavefunction, (ii) compressed
representations of the many-body wavefunction, or (iii) evolution of an
effective (downfolded) representation through Green’s functions. In this work,
we select representative cases of each of the methods and address how these
complementary approaches capture the dynamics driven by intense field
perturbations to non-equilibrium states. Under strong driving, the systems are
characterized by strong entanglement of the single particle density matrix and
natural populations approaching those of a strongly interacting equilibrium
system. We generate a representative set of results that are numerically exact
and form a basis for critical comparison of the distinct families of methods.
We demonstrate that the compressed formulation based on similarity-transformed
Hamiltonians (coupled cluster approach) is practically exact in weak fields
and, hence, weakly or moderately correlated systems. Coupled cluster, however,
struggles for strong driving fields, under which the system exhibits strongly
correlated behavior, as measured by the von Neumann entropy of the single
particle density matrix. The dynamics predicted by Green’s functions in the
(widely popular) $GW$ approximation are less accurate but improve
significantly upon the mean-field results in the strongly driven regime.

††preprint: APS/123-QED

##  I Introduction

The properties of non-equilibrium quantum systems have drawn great attention
in recent years. Driving a system can alter its properties, leading to phase
changes[1, 2], exotic states of matter[3, 4] and allowing for tuning of
material properties through continuous driving[5, 6]. Motivated by new
experimental techniques and observations, the theoretical study of these non-
equilibrium systems is also attracting wide attention[7, 8, 9, 10, 11, 12,
13]. This has lead to a wide range of theoretical techniques, that have found
success for equilibrium problems, to be extended to the non-equilibrium
regime. The available ab-initio techniques that are systematically improvable
either employ wave function or many-body Green’s function[14] as the basis of
their solution. The former are typically used in quantum chemistry, while the
latter is a suitable framework for condensed matter physics problems. There is
however lack of critical assessment of their most appropriate domain
applicability across various non-equilibrium regimes. While these approaches
are, to a large extent, complementary formulations of the time dependent
problem, they practically tackle the many-body description differently and
hence face distinct challenges which are explored in this paper.

Among wave function methods, the time dependent exact diagonalization (TD-ED)
is equivalent to the time evolution of the Schrödinger equation and hence it
provides benchmark results. Due to its cost, it is applicable to only small
finite systems. In the quantum chemistry community, it is referred to as time
dependent full configuration interaction (TD-FCI)111Note, that in ED the
system Hamiltonian is diagonalized providing all the eigenvectors, while in
TD-FCI customarily only few eigenvectors are obtained.. The time dependent
methods based on an approximate (truncated) CI expansion enjoy much success,
for example, time-dependent configuration interaction singles (TD-CIS) [16]
method is very well suited for ionization dynamics. However, such methods have
very limited applications, and active space-based CI methods are often
preferred. One of the methods in that direction is the extension of multi-
configurational time dependent Hartree approaches (MCTDH) [17] to atomic
problems such as ionization of helium [18]. Several other methods in the same
vein are also known, such as the time-dependent restricted active space CI
(TD-RASCI) [19], and the complete active space self-consistent field (TD-
CASSCF) [20]. The computational cost of those techniques scales exponentially
with the size of the Hilbert space, each particular formulation has a
different scaling prefactor, but are generally limited to small systems.

The reduced (non-exponential) scaling wavefunction-based approaches leverage
some form of compressed representation, e.g., as time dependent density matrix
renormalization group method (TD-DMRG) [21, 22, 23, 24] or time dependent
formulation of coupled cluster [25, 26] (TD-CC) methods. For these approaches,
the price of the reduction in computational scaling comes with the
establishment of the domain of most suitable applicability. TD-DMRG is usually
most successful in 1D systems, while TD-CC, following the ground-state
version, is typically suitable for a wide variety of weakly correlated
problems. [27, 28, 29, 30] One of the earliest formulations of TD-CC, of
particular interest to this work, was proposed by Arponen [7]. In recent
years, new developments in theory and implementation have emerged. [31, 32,
33, 34]

In contrast to the wavefunction methods, the methods applied to condensed
matter systems are most conveniently based upon the Green’s function (GF)
formulation of the many-body perturbation theory (MBPT)[35, 36, 37]. This
approach recasts the many-body problem onto a set of effective correlators,
most typically the one-quasiparticle GF[38], capturing the probability
amplitude associated with the propagation of a quasiparticle. As an effective
single particle quantity, GF MBPT can be efficiently implemented[39, 40]
treating systems with thousands of electrons [41, 42, 43]. The complexity is
accounted for via the space-time nonlocal self-energy that is systematically
built by expansion in the system fluctuations[44, 45, 46], making GF MBPT one
of the most powerful methods for studying electronic phenomena in large scale
systems. It is formally straightforward to extend GF MBPT to non-equilibrium
settings via the set of integrodifferential Kadanoff-Baym equations (KBEs)[47,
48, 49]. The cost of solving KBEs scales cubically with the number of
timesteps[50]. Multiple cost reduction schemes have been proposed recently[51,
52, 53, 54, 55, 56, 57, 58]. Yet, the most widely used of these is the
Hartree-Fock generalized Kadanoff-Baym ansatz (HF-GKBA)[55] which can be
formulated in linear scaling fashion[56, 59, 60, 61, 62]. In some cases the
HF-GKBA can even improve on elements of the KBEs, specifically the tendency of
the KBEs to be overdamped and reach artificial steady states[63, 64].

In this paper, we critically compare the distinct formulations of the time
dynamics and illustrate their performance on a series of numerical benchmarks
for a generalized Hubbard chain driven out of equilibrium by an external
pulse. Using TD-CI, TD-CC, the KBEs,the HF-GKBA and time dependent Hartree-
Fock(TD-HF), we study the dynamics of the density matrix and study the effects
of different strengths of non-equilibrium perturbation. In general, the TD-CI
results serve as near-exact baseline against which other methods are compared.

The paper first reviews the essential parts of the distinct formulations and
approximations in the theory section. The results then summarize the findings
for lattice problems perturbed by pulses of increasing strengths that drive
the system towards the regime which is shown to exhibit characteristics of
strong correlations. We map the capacity of the distinct formalisms to capture
the dynamics for those regimes and also for systems with short- and long-range
interactions, with and without non-local memory effects (in the context of GF
formalisms), and confirm our findings by testing systems of various sizes. The
details of the theory and implementation of the methods are in the
supplementary information[65], while the main text focuses on the analysis of
the results, which are summarized in the conclusion section.

##  II Theory

In this work, we explore representative formulations of three conceptually
distinct methodologies, configuration interactions (based on truncation of the
wavefunction expansion in the space of configurations represented by
determinants), coupled cluster approach (employing the compression of the
many-body wavefunction), and many-body Green’s function techniques (based on a
downfolded effective representation of excitation evolution). In the following
subsections we present the key underlying ideas and we leave much of the
technical details of each approach to the supplemental information[65].

###  II.1 Time-dependent configuration interaction

The time dependent full configuration interaction (TD-FCI) approach is taken
as a reference point for assessing the performance of the other methods. In
TD-FCI, we express the time-dependent Schrödinger equation in configuration
space:

$i\frac{\partial C_{I}}{\partial t}=\sum_{J}H_{IJ}(t)C_{J}$ (1)

where, $I$ and $J$ stand for all electronic configurations generated for a
given problem, and $C_{I}$ and $C_{J}$ are the coefficients of those
configurations. To solve this equation for a general time-dependent
Hamiltonian we numerically implement the time-ordering and exponential of a
large matrix, $H_{IJ}$, based on the Lanczos technique [66]. With this
approach, we construct a propagator in a small subspace based on a reference
form of the wavefunction at time $t$, i.e., $|\Psi(t)\rangle$. For example, at
$t=t_{0}$, we build a propagator based on the ground state wave function of
the time-independent Hamiltonian. Such a propagator is however valid only for
a short time. Through a regular renewed search of the “active” subspace we can
ultimately propagate the TD-FCI expansion for long times. The procedure is
numerically stable and yields unitary dynamics; the details of the
implementation are provided in the supplemental information[65].

###  II.2 Time-dependent coupled cluster

Coupled cluster (CC) is another formalism expanding the wavefunction as a
linear combination of configurations. In CC, one introduces a non-Hermitian
approximation to expectation values
$\braket{A}_{CC}=\braket{\Psi_{L}}{A}{\Psi_{R}}$, where the bra and ket
wavefunctions are defined as

$\displaystyle\ket{\Psi_{R}}$ $\displaystyle=e^{T}\ket{\Phi},$ (2a)
$\displaystyle\bra{\Psi_{L}}$
$\displaystyle=\bra{\Phi}\left(1+Z\right)e^{-T},$ (2b)

where $T=\sum_{\mu}\tau_{\mu}\gamma_{\mu}^{\dagger}$ and
$Z=\sum_{\mu}z_{\mu}\gamma_{\mu}$ are composed of particle-hole excitation
$\gamma_{\mu}^{\dagger}$ and de-excitation $\gamma_{\mu}$ operators, defined
with respect to the Slater determinant reference $\ket{\Phi}$.

In time-dependent CC, the amplitudes $\tau_{\mu}$ and $z_{\mu}$ carry all the
time-dependence, and are found by making the action integral,

$\begin{split}S&=\int dt\mathcal{L}(t)\\\ &=\int
dt\braket{\Phi}{\left(1+Z(t)\right)e^{-T(t)}\left(i\frac{\partial}{\partial
t}-H(t)\right)e^{T(t)}}{\Phi},\end{split}$ (3)

stationary with respect to variations in $\tau_{\mu}(t)$ and $z_{\mu}(t)$.
When no truncation is imposed on the operators $T$ and $Z$, CC is equivalent
to exact diagonalization, i.e., equivalent to FCI (discussed above) where the
wavefunction is constructed as a linear combination of all possible Slater
determinants within the given spin-orbital basis. In practice, however,
truncation is necessary. Fortunately, in weak-to-moderately correlated
systems, it is sufficient to include only single and double particle-hole
excitation (de-excitation) operators in $T$ ($Z$). This leads to the well
known CCSD approximation which is used for all the results and discussion
presented in this paper.

It is noteworthy that due to the exponential parameterization, the ket
wavefunction in truncated CC includes contributions from all Slater
determinants, although these contributions are factorized in terms of the
lower-rank parameters. This is in contrast with truncated CI, where the cut-
off is imposed explicitly on the number of determinants. At the same time, the
asymmetric nature of CC expectation values leads to a violation of the
variational principle. Consequently, in TD-CC, real observable quantities may
develop un-physical or complex parts. However, the unphysical terms often
serve as a diagnostic tool for the CCSD approximation. In the moderately
interacting regimes, where CC works well, the expectation values are well
behaved. The increasing magnitude of unphysical expectation values generally
coincides with the onset of strong correlations, where CC is known to fail.
This point will be discussed in the sections below. Further details regarding
the derivation and implementation of the TD-CC equations are provided in the
supplemental information. [65]

###  II.3 Green’s function methods

Finally we employ two approaches based on the time evolution of the one-
particle Green’s function (GF), i.e., the time dependent probability amplitude
associated with propagation of a quasiparticle. At non-equilibrium, the GF is
an explicit function of two time arguments, and their respective ordering on
the real time axis (as well as imaginary time axis for the finite temperature
formalisms). The first are the Kadanoff-Baym equations, representing a set of
theoretically exact integro-differential equations for the two-time GF, and
are given by[67]

$\begin{split}[-\partial_{\tau}-h]G^{\mathrm{M}}(\tau)&=\delta(\tau)+\int_{0}^{\beta}d\bar{\tau}\Sigma^{\mathrm{M}}(\tau-\bar{\tau})G^{\mathrm{M}}(\bar{\tau}),\\\
i\partial_{t_{1}}G^{\rceil}(t_{1},-i\tau)&=h^{\textrm{HF}}(t_{1})G^{\rceil}(t_{1},-i\tau)+I^{\rceil}(t_{1},-i\tau),\\\
-i\partial_{t_{2}}G^{\lceil}(-i\tau,t_{2})&=G^{\lceil}(-i\tau,t_{2})h^{\textrm{HF}}(t_{2})+I^{\lceil}(-i\tau,t_{1}),\\\
i\partial_{t_{1}}G^{\lessgtr}(t_{1},t_{2})&=h^{\textrm{HF}}(t)G^{\lessgtr}(t_{1},t_{2})+I_{1}^{\lessgtr}(t_{1},t_{2}),\\\
-i\partial_{t_{2}}G^{\lessgtr}(t_{1},t_{2})&=G^{\lessgtr}(t_{1},t_{2})h^{\textrm{HF}}(t_{2})+I_{2}^{\lessgtr}(t_{1},t_{2}),\\\
\end{split}$ (4)

Here $\Sigma$ is a space-time nonlocal effective potential embodying all many-
body interactions; $G^{\mathrm{M},<,>,\rceil,\lceil}$ are various components
of the Green’s function that depend on which of the time arguments are in real
time versus in imaginary time and the $I^{<,>,\rceil,\lceil}$ are integrals
that account for many-body correlation effects in the system. A full
description and discussion of these quantities is given in the supplemental
information[65].

The KBEs are approximate only through the choice of the self-energy. In this
work, we focus on the $GW$ approximation, representing the arguably most
popular choice for simulations of condensed matter systems. In essence, this
approach accounts for dynamical density-density interactions, which dominate
the correlation effects in weakly and moderately correlated systems (such as
semiconductors[44, 68, 14]. The combination of integrals in the equations
above and their two-time nature means solving the KBEs scales cubically in the
number of time steps, making these calculations prohibitively expensive.
Furthermore, the KBEs often suffer from issues of artificial damping coming
from the approximate self-energy[64, 63].

Given the cost of KBEs propagation, they often solved only approximately, and
the most common route is to employ the Hartree-Fock Generalized Kadanoff-Baym
ansatz (HF-GKBA). It assumes a particular self-energy for the time diagonal
component of the GF and takes the off-diagonal self-energy to be approximated
by only the bare exchange interactions (hence the Hartree-Fock (HF)
approximation in its name). By using the HF form of the time off-diagonals, it
is possible to derive an ordinary differential equation (ODE) scheme for the
evolution of the equal-time component of the GF. This translates to a time-
linear formalism that recasts the equation of motion using the explicit form
of the two particle GF instead of making use of $\Sigma$ [69, 56]. The
equations of motion for this scheme are given by:

$\begin{split}i\partial_{t}G^{<}_{ij}(t)&=[h^{\textrm{HF}}(t),G^{<}(t)]_{ij}+[I+I^{\dagger}]_{ij}(t),\\\
I_{ij}(t)&=-i\sum_{klp}w_{iklp}(t)\mathcal{G}_{lpjk}(t),\\\
i\partial_{t}\mathcal{G}_{ijkl}(t)&=[h^{(2),\textrm{HF}}(t),\mathcal{G}(t)]_{ijkl}\\\
&\hskip 56.9055pt+\Psi_{ijkl}(t)+\Pi_{ijkl}-\Pi_{lkji}^{*}.\end{split}$ (5)

Here $\mathcal{G}(t)$ is a two particle propagator, $\Psi_{ijkl}$ accounts for
pair correlations built up due to two-particle scattering events and
$\Pi_{ijkl}$ accounts for polarization effects[69]. While approximate, the HF-
GKBA does not suffer from artificial damping as the KBEs do and thus in some
cases offers improved results over the KBEs at a much reduced cost. This
further enables the use of more advanced forms of $\Sigma$, the study of
larger systems, and the ability to perform longer time evolution. A more
detailed description of the HF-GKBA and the ODE scheme are given in the
supplemental information[65].

###  II.4 Model Systems

To benchmark different time-dependent methods, we use a generalization of the
one-dimensional Hubbard model, which is driven out of equilibrium with the
help of an electric-field pulse. The Hamiltonian for the system is written as:

$\begin{split}\mathcal{H}&=-J\sum_{\langle
i,j\rangle\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\\\
&\hskip
71.13188pt+V\sum_{i}(-1)^{i}n_{i}+\sum_{ij}h^{\mathrm{N.E}}_{ij}c_{i}^{\dagger}c_{j}.\end{split}$
(6)

The first two terms are the usual nearest neighbor hopping and onsite
interaction of the Hubbard model. The third term turns the system into a
bipartite lattice and serves to open up a trivial gap in the system. This is
necessary since the HF-GKBA prepares the initial state using adiabatic
switching. The adiabatic theorem requires a gapped system to correctly prepare
the HF-GKBA initial state. Such gapped models have been used in previous
studies of the HF-GKBA[57, 58]. We choose $J=1$ and express all interaction
strengths in units of $J$. Further, all the results and discussion assume half
filling, i.e., one electron per site.

For each method the system is prepared in the corresponding ground state
before being excited with a time dependent field modeled under the dipole
approximation and given by

$h^{\mathrm{N.E}}_{ij}(t)=\delta_{ij}\left(\frac{N-1}{2}-i\right)E\mathrm{e}^{-\frac{(t-t_{0})^{2}}{2\sigma^{2}}}.$
(7)

Here $N$ is the number of lattice sites in the model, $E$ is the strength of
the field, $\sigma$ determines the temporal width of the pulse and $t_{0}$ is
the pulse midpoint.

##  III Results

###  III.1 Effect of excitation strength

First we investigate the effect of the strength of the perturbation by varying
the parameter $E$ in our model. This will capture the behavior of each method
as we go from linear response type excitations up to strongly perturbed
regimes. In Fig. 1 we show the dynamics of the electronic dipole for the
system in equation (6) and (7) with $E=1.0J$, $2.0J$, $3.0J$ and $5.0J$, as
computed with each of our methods. For this section, we show results obtained
for the number of sites $N=12$ and we choose an interaction strength of
$U=1.0J$. This puts the kinetic and Coulomb energy scales at the same level.
Here each method captures the $U=1.0$ ground states very well. Thus, this
setup allows us to see how each method performs in non-equilibrium without
having to worry about the effects of a poorly captured ground state. While we
have explored various parameters, the findings below are representative of the
behaviors observed. Additional results are discussed in the supplemental
information[65].

We will start this section by discussing the results produced by the reference
TD-CI approach, which is numerically exact. These results will serve as our
benchmark to compare against other methods. For $E=1.0J$, we see two dominant
types of oscillation: a strong low-frequency contribution that manifests as
long wavelength oscillations in the dipole and a weak high-frequency
contribution. This behavior continues with increasing $E,$ however, for strong
perturbation strengths, we see that the magnitude of both types of
oscillations diminishes, and a yet lower frequency mode begins to dominate.
This continues until the dynamics become almost flat when we come to the
bottom panel of Fig. 1 with $E=5.0J$. Further results and analysis to
understand these exact trends in dipole moment are discussed in Sec. III.2 and
IV.

We continue this section by comparing the results produced by our other
methods.

For completeness, we include the TD-HF trajectory, which serves as a baseline
for the types of expansion techniques discussed below (either CC or those
based on the GFs, with HF being their underlying reference Hamiltonian around
which the self-energy is expanded). TD-HF is based solely on a time evolution
using a single determinant. For $E=1.0J$ we see TD-HF performs well, capturing
the amplitude and frequency of the oscillations of the TD-CI. It suffers from
a similar issue the HF-GKBA has in that it has a slight offset in frequency
that causes it to move out of phase with the exact result as the time
evolution progresses. Interestingly for $E=2.0J$ the TD-HF actually improves
upon the results given by HF-GKBA and the KBEs. This fortuitous improvement is
likely due to the self-energy over-damping the MBPT results. As $E$ is
increased again we see the result worsen. For $E=3.0J$ TD-HF is in agreement
up to $t\approx 10J^{-1}$ and for $E=5.0J$ to around $t=8J^{-1}$. After these
times however the TD-HF result differs drastically, although still oscillates
around the correct result. Judging by the results, it seems now the TD-HF
result is not damped enough, due to the lack of a self-energy and thus
oscillates with much greater amplitude than the benchmark results. In the next
paragraphs, we discuss the CC and GF-based techniques, which can be thought of
as expansions around this “correlation-free” TD-HF dynamics.

The TD-CC approach is in excellent agreement with TD-CI for the lowest two $E$
values ($E=1.0J,2.0J$). Note that this is despite the relatively complex
dynamics. For increased strength (starting with $E=3.0J$), TD-CC captures the
dynamics very closely only up until $t\approx 10J^{-1}$, but deviates from the
exact solution soon after. This illustrates a growing error in TD-CC for long-
time evolutions at large excitation strengths despite the initial TD-CC
dynamics for $E=3.0J$ being in reasonable qualitative agreement with TD-CI.
Note that without the TD-CI reference, it is not immediately clear from Fig. 1
that TD-CC dynamics is erroneous. For $E=4.0J,5.0J$, soon after the system is
driven out of equilibrium (around $t\approx 5J^{-1}$), the CC approach starts
to break down, and the results become unphysical as the occupation numbers
develop a large imaginary contribution. In fact, the disagreement of CC with
TD-CI in the presence of strong perturbation can be anticipated by looking at
Fig. S6 of SI, [65] where we plot the imaginary part of CC estimates for the
occupation number on the first site in the Hubbard lattice. It is evident that
even for $E=3.0J$ at long time scales, and especially for $E=4.0J,\,5.0J$, TD-
CC leads to large unphysicalities, signaling its breakdown. At the same time,
keeping track of imaginary parts or inconsistent behavior in expectation
values in TD-CC can serve as a diagnostic tool to measure the efficacy of the
coupled cluster approximation itself.

Next, we turn to the MBPT methods, starting with the solution for the KBEs
using the $GW$ self-energy. For $E=1.0J$ the agreement is very strong up to
$t\approx 15J^{-1}$. After this, while TD-CC improves in the quantitative
results, the KBEs still provide qualitatively correct results. They capture
the low frequency mode of oscillation extremely well, however the high
frequency part is damped out almost entirely by $t\approx 25J^{-1}$. Moving to
$E=2.0J$, the time up to which we see perfect agreement is shortened to
$t\approx 10J^{-1}$. Additionally, unlike the $E=1.0J$ case where reasonable
qualitative agreement held until $t=50J^{-1}$, here we see that after
$t\approx 20J^{-1}$ the KBE results differ significantly from those produced
by TD-CI and TD-CC. The damping in the KBEs is more pronounced now, and the
dynamics are close to flat. This damping represents a well-known issue with
the KBEs and will be discussed later in section IV. This trend continues as we
continue to increase the excitation strength. For $E=3.0J$, perfect agreement
with TD-CI now only holds to $t\approx 10J^{-1}$, around the same as the TD-CC
result. The KBE damping kicks in almost immediately after this at $t\approx
12J^{-1}$. After this the KBE dynamics are essentially flat. Interestingly,
the exact TD-CI result seems to oscillate around the stationary value that the
KBEs reach. For $E=5.0J$ the KBE results fail almost immediately by strongly
suppressing any oscillations at or after $t\approx 5J^{-1}$. Again, we see the
result is damped (leading to only a slowly time-varying dipole moment which
monotonically increases within the simulation window). However, compared to
the TD-CI result, this is not qualitatively wrong despite the magnitude of the
dipole is now much further from the TD-CI result compared to the results in
Fig. 1(a)-(c).

We now turn to the HF-GKBA, which truncates certain memory effects compared to
the KBEs. Similar to KBEs, we see an excellent agreement with TD-CI for
$E=1.0J$ up until $t\approx 15J^{-1}$. After this, the HF-GKBA still shows a
strong qualitative and quantitative agreement compared to the benchmark
results, contrasting with the overdamped behavior of KBEs. The HF-GKBA
captures the low frequency mode of oscillation very well. For longer times,
the high frequency oscillations from the TD-CI and TD-CC results are offset.
Additionally, the magnitude of these high frequency oscillations is diminished
relative to the benchmark results. However, the HF-GKBA does systematically
improve on the damping seen in the KBE result. For $E=2.0J$ the HF-GKBA
follows a very similar trend to the KBEs. The main difference is, again, a
reduction in the damping seen in the KBE dynamics. Both follow roughly the
same trajectory, which appears to be centered around the main low frequency
oscillations of TD-CI and TD-CC. Moving to $E=3.0J$ the HF-GKBA result now has
a noticeable offset from the KBE solution, but keeps the same trend of
improving upon the overdamping. After $t\approx 10J^{-1}$ the HF-GKBA no
longer captures the benchmark results qualitatively or quantitatively.
Furthermore, the HF-GKBA oscillates minorly around a single value (somewhat
similarly to the KBEs). However, the value the HF-GKBA result seems to center
around is lower than the central values of the TD-CI, TD-CC, and the KBEs
which are more or less equivalent. Turning finally to Fig. 1(d) with
$E=5.0J^{-1}$ we see the HF-GKBA actually captures the qualitative behavior of
the TD-CI dipole well. The HF-GKBA displays a slightly larger magnitude of
oscillation than those of TD-CI, however the mean values are very close and it
offers a markedly significant improvement over the TD-CC and KBE results.

Figure 1: Dynamics of the electronic dipole system described by Eq. (6) and
(7) for various values of $E$ parameter. a) $E=1.0J$, b) $2.0J$, c) $3.0J$ d)
$5.0J$.

###  III.2 Exciting to strongly correlated non-equilibrium regime

In this section, we further investigate the results discussed in section
III.1. As the TD-CI results provide access to the (nearly exact) wavefunction
trajectory, we use it to analyze the time evolution in greater depth. In the
next section (IV), we rely on this analysis to illustrate the origins of
failure in each method in capturing the proper physics of non-equilibrium
systems.

In Fig. 2(a)-(c), we show the natural occupations of the single particle
density matrix for the model discussed in section II.4 and section III.1 for
$E=1.0J$, $E=2.0J$ and $E=5.0J$. The natural occupations correspond to the
eigenvalues of time dependent single particle density matrix expressed in
spin-orbitals, taken from TD-CI. Our analysis is motivated by equilibrium,
ground-state physics of the underlying half-filled Hubbard model, where
natural occupations approach $\lambda_{i}\rightarrow 0.5$ as the system
becomes strongly correlated in the $U\rightarrow\infty$ limit (see Fig. S3 in
supplemental information[65]). By examining natural occupations, one can
therefore assess how correlated the system is. Here we propose this quantity
as a measure of correlation in the non-equilibrium setting.

In Fig. 2 (d), we also look at the von-Neumann entropy for a scan of $E$
values from $1.0J$ to $5.0J$. The entropy is computed as,

$\mathcal{S}[\rho]=\frac{2}{N}\mathrm{Tr}\left[\rho\log_{2}(\rho)\right].$ (8)

Note that for the ground-state, as $U\to\infty$, the entropy approaches unity.

For $E=1.0J$, in Fig. 2(a), we see very little deviation from the natural
occupations. This corresponds to near zero entropy of the single particle
density matrix. Each of our methods matches the benchmark perfectly in
equilibrium for the results in Fig 1. This demonstrates that when the natural
occupations remain close to the equilibrium values the four methods capture
well the non-equilibrium TD-CI results qualitatively and even quantitatively.

As the excitation strength is increased to $E=2.0J$, we see a more noticeable
difference from the equilibrium result. Though the oscillations of the natural
occupations are still relatively small. The entropy reflects this more clearly
where it is only around $\mathcal{S}[\rho]=0.1$, indicating weak mixing
between the natural orbitals. For the dynamics in Fig. 1, the results produced
by TD-CC are still in perfect agreement with the TD-CI result. For the same
conditions, the MBPT methods appear over damped but they still capture
elements of the TD-CI result.

For $E=5.0J$ we see a strong mixing of the natural occupations, nearing the
$\lambda_{i}[\rho]=0.5$ limit, corresponding to an entropy of
$\mathcal{\mathcal{S}[\rho]}\approx 0.9$. In an equilibrium picture this would
correspond to a strongly correlated system, with high $U$. Such a transition
was observed for all system sizes and interaction range types. Consequently,
we consider this to be a representative case of a system driven by a strong
external perturbation from a “weakly correlated” ground state to a “strongly
correlated” excited state.

Understandably, such a physical regime poses difficulty for methods suited for
weakly correlated systems, such as CC or the MBPT expansion. Indeed, TD-CC
results fail completely in this regime, as do the KBE results. In contrast the
HF-GKBA offers an improvement in this regime and does a good job at capturing
qualitative and quantitative of the TD-CI results. In Fig. S5 of the
supplemental information we show the entropy for different $E$ for
$N_{s}=8,12$ and $16$ and see similar behavior for each system size[65].

Figure 2: Left: Natural occupations of the single particle density matrix
taken from TD-CI for the model discussed in section II.4 and section III.1 for
a) $E=1.0J$, b) $E=2.0J$ and c) $E=5.0J$. Right: The entropy for a scan of $E$
values from $1.0J$ to $5.0J$

##  IV Discussion and conclusions

We now return to discuss the implications of the results we have described in
section III.1 and III.2. Each of the methods we test perform reasonably well
in weakly excited systems. For pulse strengths up to $E=2.0J$ TD-CC gives
perfect agreement compared to TD-CI. However, as stronger perturbations are
applied to drive the system out of equilibrium, TD-CC breaks down. This is
because CC wavefunction is expanded around the time-independent, ground-state
Hartree-Fock reference which becomes increasingly bad as the system evolves
after a strong time-dependent perturbation.

Although, even for $E=3.0J$ TD-CC offers a good improvement over the HF-GKBA
and the KBEs. Interestingly, at $E=2.0J$ we observe that TD-HF also improves
upon the results for the time-dependent dipole produced by HF-GKBA and the
KBEs. Note that this implies that the trajectories of equal-time observables
depending on the density matrix (such as the dipole studied here) are
qualitatively well captured even by a mean-field technique. Further, this
leads to the interesting question of whether in dynamical systems it is always
better to include self-energy effects, or whether there are regimes in which
the inclusion of the (non-exact) self-energy can have detrimental effects by
overdamping the system dynamics. This is a question that deserves further
investigation with additional forms of the self-energy beyond merely the
density-density response included in the $GW$ formalism, but also including
higher order fluctuations[45, 46].

In the $E=1.0J$ case the Kadanoff-Baym equations qualitatively capture the
dynamics after excitation. However, as is commonly seen[64, 63] they suffer
from severe overdamping that removes the high frequency component of the
dipole almost immediately. For $E=1.0J$ the HF-GKBA also captures the
qualitative features of the dipole oscillations and even offers improvement to
the KBEs by removing some of the damping. As we increase $E$ further, the full
KBEs continue to be overdamped, while the HF-GKBA improves the result, adding
back in some of the high frequency oscillation modes lost in the KBEs.

Interestingly, as we show in Fig. S4 of the supplemental information[65], when
long range interactions are included in the Hamiltonian the result of the HF-
GKBA are improved drastically. This improvement is likely related to the
suitability of a self-energy approximation for a given problem. The $GW$
approximation, used here, assumes the dominant contribution to the self-energy
is from the screened coulomb interaction, or equivalently from the induced
charge density fluctuations. Due to the extremely localized nature of the
Hubbard model screening is relatively weak, meaning a different self-energy
approximation may be more appropriate. A long range interacting model
naturally has more potential for screening, thus making the $GW$ approximation
a more relevant contribution to the self-energy, and so improving the
agreement with TD-CI.

This observation may also help to explain the improvement of the HF-GKBA for
$E=5.0J$. In this highly excited regime, the particles are much more
delocalized. This increases the screening and we propose that this lead to a
regime in which the physics accounted for by the $GW$ approximation becomes
more dominant, leading to reasonably well captured dynamics. Indeed, $GW$ can
be interpreted as dynamically screened exchange interaction, i.e., its
difference compared to TD-HF stems purely from adding correlations effects due
to the density fluctuations which likely dominate in this regime, over other
types of correlations. It is in this regime that CC shows its most significant
failing, showing results that are drastically different from the CI result
both quantitatively and qualitatively. Here, the main difference between GF
MBPT and CC stems from resummations of the screened interactions, which are
not included in CC. TD-CC can possibly be improved via orbital optimization in
the pursuit of renormalizing the interactions, as well as by including
additional higher-order excitation terms (beyond the CCSD level).[31, 32, 70,
33]

From these results, it is clear that the TD-CC performs excellently for
relatively weak perturbations and improves upon MBPT, especially in the finer
details. However, in strongly perturbed regimes, TD-CC breaks down and no
longer provides results that match TD-CI, even at the qualitative level.
Notably, for regimes in which the system is strongly perturbed, the HF-GKBA
offers results that are very close to those produced by TD-CI, and are much
more numerically stable than those produced by TD-CC. This points to the
possibility that the HF-GKBA and TD-CC have somewhat complimentary regimes of
applicability, however definitive statements about this observation will
require further exploration of various self-energy formulations. As is often
seen the KBEs show no significant improvement upon the HF-GKBA, furthermore in
this particular example due to the severe overdamping the HF-GKBA outperforms
the KBEs. This points to the fact that inclusion of the time-nonlocal
components of the self-energy, which dominate the “quantum memory” of KBEs, is
practically detrimental to the quality of the predicted dynamics, despite its
formal correctness. This question warrants further exploration, especially in
connection to applications of more involved self-energy formulations, beyond
$GW$.

We link the observed performance of the CC and GF methods (within the given
approximations) to the effective strength of correlations induced by the non-
equilibrium. In section III.2, we looked at the different regimes brought
about by varying the pulse strength. The entropy of the density matrix, given
in equation (8), is here used to quantify correlations in the system[71]. In
Fig. 2 we showed that for stronger pulses the natural populations become
closer to $\lambda=0.5$, which is consistent with those of a strongly
interacting equilibrium system as shown in Fig. S3 of the supplemental
information for the model in equation (6) with $N_{s}=6$[65]. Further, in in
Fig. 1 (d) the dynamics of the dipole are very slowly varying in time. This is
consistent with the type of dynamics one would see in a very strongly
interacting system due suppression of the kinetic energy term. We interpret
this regime as non-equilibrium induced correlations.

In this work we have only investigated properties related to the density
matrix, namely the dipole and the von Neumann entropy. Capturing the dynamics
of the density matrix is an important benchmark for any time dependent method,
and being able to predict these quantities is valuable for quantitative
simulation of carrier dynamics in materials or charge dynamics in quantum
chemistry systems. On the other hand, other measurable quantities of great
importance for non-equilibrium quantum physics, such as the the time-dependent
spectral function, require knowledge of the time-non-local terms. We are
actively investigating the ability of each of these methods to capture the
time dependent spectral properties of non-equilibrium systems, as well as
looking into efficient ways of improving existing approaches[72].

A further extension of this study, also underway, comes in terms of the system
sizes studied. Despite the huge reduction in cost afforded in TD-CI compared
to full diagonalization, we are still limited in the system size due to the
still large size of the constructed subspace. One route we are actively
pursuing is implementing a time dependent adaptive sampling CI approach(TD-
ASCI). This approach improves upon traditional CI by adding a step that
weights wavefunctions making up the subspace, based on estimating how much of
a contribution they will make to the subspace. This can vastly reduce the
subspace size thus opening the road to even larger systems. This presents a
route to provide accurate benchmark data in large scale or multi-band periodic
systems where such data is extremely difficult to obtain.

In conclusion, the numerical benchmarks presented in this work help to
practically assess the performance of variety of time-dependent numerical
methods under different strengths of perturbation. We demonstrate that TD-CC
is an excellent theory for moderate coupling and weak to moderately strong
perturbations, but that it becomes unreliable for strong perturbations. We
have further seen the KBEs fail to produce results that improve upon TD-CC.
For weak excitations the HF-GKBA produces worse results than TD-CC but for
strong excitations the HF-GKBA gives results that are qualitatively and
quantitatively very close to the true result. Thus, we surmise that a
multilayered methodology employing the complementary strengths of TD-CC and
HF-GKBA will likely provide accurate results in a wide variety of non-
equilibrium regimes.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research and Office
of Basic Energy Sciences, Scientific Discovery through Advanced Computing
(SciDAC) program under Award Number DE-SC0022198. This research used resources
of the National Energy Research Scientific Computing Center, a DOE Office of
Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award
BES-ERCAP0029462

## References

  * Beebe  _et al._ [2017] M. R. Beebe, J. M. Klopf, Y. Wang, S. Kittiwatanakul, J. Lu, S. A. Wolf, and R. A. Lukaszew, Time-resolved light-induced insulator-metal transition in niobium dioxide and vanadium dioxide thin films, Opt. Mater. Express 7, 213 (2017).
  * Disa  _et al._ [2023] A. S. Disa, J. Curtis, M. Fechner, A. Liu, A. von Hoegen, M. Först, T. F. Nova, P. Narang, A. Maljuk, A. V. Boris, B. Keimer, and A. Cavalleri, Photo-induced high-temperature ferromagnetism in YTiO3, Nature 617, 73 (2023).
  * Dong  _et al._ [2021] S. Dong, M. Puppin, T. Pincelli, S. Beaulieu, D. Christiansen, H. Hübener, C. W. Nicholson, R. P. Xian, M. Dendzik, Y. Deng, Y. W. Windsor, M. Selig, E. Malic, A. Rubio, A. Knorr, M. Wolf, L. Rettig, and R. Ernstorfer, Direct measurement of key exciton properties: Energy, dynamics, and spatial distribution of the wave function, Natural Sciences 1, e10010 (2021).
  * Bao  _et al._ [2022] C. Bao, P. Tang, D. Sun, and S. Zhou, Light-induced emergent phenomena in 2D materials and topological materials, Nature Reviews Physics 4, 33 (2022).
  * Nuske  _et al._ [2020] M. Nuske, L. Broers, B. Schulte, G. Jotzu, S. A. Sato, A. Cavalleri, A. Rubio, J. W. McIver, and L. Mathey, Floquet dynamics in light-driven solids, Phys. Rev. Res. 2, 043408 (2020).
  * Zhou  _et al._ [2023] S. Zhou, C. Bao, B. Fan, F. Wang, H. Zhong, H. Zhang, P. Tang, W. Duan, and S. Zhou, Floquet engineering of black phosphorus upon below-gap pumping, Phys. Rev. Lett. 131, 116401 (2023).
  * Arponen [1983] J. Arponen, Variational principles and linked-cluster exp S expansions for static and dynamic many-body problems, Annals of Physics 151, 311 (1983).
  * Ishikawa and Sato [2015] K. L. Ishikawa and T. Sato, A Review on Ab Initio Approaches for Multielectron Dynamics, IEEE Journal of Selected Topics in Quantum Electronics 21, 1 (2015).
  * Li  _et al._ [2020] X. Li, N. Govind, C. Isborn, A. E. DePrince, and K. Lopata, Real-Time Time-Dependent Electronic Structure Theory, Chem. Rev. 120, 9951 (2020).
  * Sun  _et al._ [2021] J. Sun, C.-W. Lee, A. Kononov, A. Schleife, and C. A. Ullrich, Real-Time Exciton Dynamics with Time-Dependent Density-Functional Theory, Phys. Rev. Lett. 127, 077401 (2021).
  * Karlsson  _et al._ [2021] D. Karlsson, R. van Leeuwen, Y. Pavlyukh, E. Perfetto, and G. Stefanucci, Fast Green’s Function Method for Ultrafast Electron-Boson Dynamics, Phys. Rev. Lett. 127, 036402 (2021).
  * Perfetto  _et al._ [2022] E. Perfetto, Y. Pavlyukh, and G. Stefanucci, Real-Time $GW$: Toward an ab initio description of the ultrafast carrier and exciton dynamics in two-dimensional materials, Phys. Rev. Lett. 128, 016801 (2022).
  * Sverdrup Ofstad  _et al._ [2023] B. Sverdrup Ofstad, E. Aurbakken, Ø. Sigmundson Schøyen, H. E. Kristiansen, S. Kvaal, and T. B. Pedersen, Time-dependent coupled-cluster theory, WIREs Computational Molecular Science 13, e1666 (2023).
  * Martin  _et al._ [2016a] R. M. Martin, L. Reining, and D. M. Ceperley,  _Interacting electrons_ (Cambridge University Press, 2016).
  * Note [1] Note, that in ED the system Hamiltonian is diagonalized providing all the eigenvectors, while in TD-FCI customarily only few eigenvectors are obtained.
  * Rohringer  _et al._ [2006] N. Rohringer, A. Gordon, and R. Santra, Configuration-interaction-based time-dependent orbital approach for ab initio treatment of electronic dynamics in a strong optical laser field, Phys. Rev. A 74, 043420 (2006).
  * Lode  _et al._ [2020] A. U. J. Lode, C. Lévêque, L. B. Madsen, A. I. Streltsov, and O. E. Alon, Colloquium: Multiconfigurational time-dependent hartree approaches for indistinguishable particles, Rev. Mod. Phys. 92, 011001 (2020).
  * Hochstuhl and Bonitz [2011] D. Hochstuhl and M. Bonitz, Two-photon ionization of helium studied with the multiconfigurational time-dependent Hartree–Fock method, The Journal of Chemical Physics 134, 084106 (2011), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3553176/14905108/084106_1_online.pdf .
  * Hochstuhl and Bonitz [2012] D. Hochstuhl and M. Bonitz, Time-dependent restricted-active-space configuration-interaction method for the photoionization of many-electron atoms, Phys. Rev. A 86, 053424 (2012).
  * Sato and Ishikawa [2013] T. Sato and K. L. Ishikawa, Time-dependent complete-active-space self-consistent-field method for multielectron dynamics in intense laser fields, Phys. Rev. A 88, 023402 (2013).
  * Daley  _et al._ [2004] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Zwolak and Vidal [2004] M. Zwolak and G. Vidal, Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm, Phys. Rev. Lett. 93, 207205 (2004).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Crawford and Schaefer [2000] T. D. Crawford and H. F. Schaefer, An Introduction to Coupled Cluster Theory for Computational Chemists, in  _Reviews in Computational Chemistry_, edited by K. B. Lipkowitz and D. B. Boyd (John Wiley & Sons, Inc., 2000) pp. 33–136.
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, Coupled-cluster theory in quantum chemistry, Rev. Mod. Phys. 79, 291 (2007).
  * Wälz  _et al._ [2012] G. Wälz, D. Kats, D. Usvyat, T. Korona, and M. Schütz, Application of Hermitian time-dependent coupled-cluster response Ansätze of second order to excitation energies and frequency-dependent dipole polarizabilities, Phys. Rev. A 86, 052519 (2012).
  * Koulias  _et al._ [2019] L. N. Koulias, D. B. Williams-Young, D. R. Nascimento, A. E. I. DePrince, and X. Li, Relativistic Real-Time Time-Dependent Equation-of-Motion Coupled-Cluster, J. Chem. Theory Comput. 15, 6617 (2019).
  * Skeidsvoll  _et al._ [2020] A. S. Skeidsvoll, A. Balbi, and H. Koch, Time-dependent coupled-cluster theory for ultrafast transient-absorption spectroscopy, Phys. Rev. A 102, 023115 (2020).
  * Skeidsvoll and Koch [2023] A. S. Skeidsvoll and H. Koch, Comparing real-time coupled-cluster methods through simulation of collective Rabi oscillations, Phys. Rev. A 108, 033116 (2023).
  * Kvaal [2012] S. Kvaal, Ab initio quantum dynamics using coupled-cluster, J. Chem. Phys. 136, 194109 (2012).
  * Sato  _et al._ [2018] T. Sato, H. Pathak, Y. Orimo, and K. L. Ishikawa, Communication: Time-dependent optimized coupled-cluster method for multielectron dynamics, J. Chem. Phys. 148, 051101 (2018).
  * Pathak  _et al._ [2020] H. Pathak, T. Sato, and K. L. Ishikawa, Time-dependent optimized coupled-cluster method for multielectron dynamics. II. A coupled electron-pair approximation, The Journal of Chemical Physics 152, 124115 (2020).
  * Peng  _et al._ [2021] R. Peng, A. F. White, H. Zhai, and G. Kin-Lic Chan, Conservation laws in coupled cluster dynamics at finite temperature, The Journal of Chemical Physics 155, 044103 (2021).
  * Reining [2018] L. Reining, The $GW$ approximation: content, successes and limitations, WIREs Computational Molecular Science 8, e1344 (2018).
  * Onida  _et al._ [2002] G. Onida, L. Reining, and A. Rubio, Electronic excitations: density-functional versus many-body green’s-function approaches, Rev. Mod. Phys. 74, 601 (2002).
  * Aulbur  _et al._ [2000] W. G. Aulbur, L. Jönsson, and J. W. Wilkins, Quasiparticle calculations in solids, Journal of Physics C: Solid State Physics 54, 1 (2000).
  * Martin  _et al._ [2016b] R. Martin, L. Reining, and D. Ceperley,  _Interacting Electrons_ (Cambridge University Press, 2016).
  * Wilhelm  _et al._ [2021] J. Wilhelm, P. Seewald, and D. Golze, Low-scaling gw with benchmark accuracy and application to phosphorene nanosheets, Journal of Chemical Theory and Computation 17, 1662 (2021), pMID: 33621085, https://doi.org/10.1021/acs.jctc.0c01282 .
  * Neuhauser  _et al._ [2014] D. Neuhauser, Y. Gao, C. Arntsen, C. Karshenas, E. Rabani, and R. Baer, Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic $GW$ Approach, Phys. Rev. Lett. 113, 076402 (2014).
  * Vlček  _et al._ [2018] V. Vlček, W. Li, R. Baer, E. Rabani, and D. Neuhauser, Swift $GW$ beyond 10,000 electrons using sparse stochastic compression, Phys. Rev. B 98, 075107 (2018).
  * Gao  _et al._ [2016] W. Gao, W. Xia, X. Gao, and P. Zhang, Speeding up gw calculations to meet the challenge of large scale quasiparticle predictions, Scientific Reports 6, 36849 (2016).
  * Govoni and Galli [2015] M. Govoni and G. Galli, Large scale gw calculations, Journal of Chemical Theory and Computation 11, 2680 (2015), pMID: 26575564, https://doi.org/10.1021/ct500958p .
  * Hedin [1965] L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem, Phys. Rev. 139, A796 (1965).
  * Vlcek [2019] V. Vlcek, Stochastic vertex corrections: Linear scaling methods for accurate quasiparticle energies, Journal of chemical theory and computation 15, 6254 (2019).
  * Mejuto-Zaera and Vlček [2022] C. Mejuto-Zaera and V. Vlček, Self-consistency in $GW\mathrm{\Gamma}$ formalism leading to quasiparticle-quasiparticle couplings, Phys. Rev. B 106, 165129 (2022).
  * Stefanucci and van Leeuwen [2013] G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, 2013).
  * Kadanoff and Baym [1961] L. Kadanoff and G. Baym,  _Quantum Statistical Mechanics_ (The Benjamin/Cummings Publishing Company, New York, 1961).
  * Schwinger [1961] J. Schwinger, Brownian Motion of a Quantum Oscillator, Journal of Mathematical Physics 2, 407 (1961), https://pubs.aip.org/aip/jmp/article-pdf/2/3/407/19046507/407_1_online.pdf .
  * Bonitz [2015] M. Bonitz,  _Quantum Kinetic Theory_ (Springer International Publishing Switzerland, 2015).
  * Kaye and Golež [2021] J. Kaye and D. Golež, Low rank compression in the numerical solution of the nonequilibrium Dyson equation, SciPost Phys. 10, 091 (2021).
  * Kaye and U. R. Strand [2023] J. Kaye and H. U. R. Strand, A fast time domain solver for the equilibrium Dyson equation, Advances in Computational Mathematics 49, 10.1007/s10444-023-10067-7 (2023).
  * Yin  _et al._ [2021] J. Yin, Y.-h. Chan, F. da Jornada, D. Qiu, C. Yang, and S. G. Louie, Analyzing and predicting non-equilibrium many-body dynamics via dynamic mode decomposition, arXiv preprint arXiv:2107.09635 (2021).
  * Stahl  _et al._ [2022] C. Stahl, N. Dasari, J. Li, A. Picano, P. Werner, and M. Eckstein, Memory truncated Kadanoff-Baym equations, Phys. Rev. B 105, 115146 (2022).
  * Lipavský  _et al._ [1986] P. Lipavský, V. Špička, and B. Velický, Generalized Kadanoff-Baym ansatz for deriving quantum transport equations, Phys. Rev. B 34, 6933 (1986).
  * Schlünzen  _et al._ [2020] N. Schlünzen, J.-P. Joost, and M. Bonitz, Achieving the Scaling Limit for Nonequilibrium Green Functions Simulations, Phys. Rev. Lett. 124, 076601 (2020).
  * Reeves  _et al._ [2023a] C. C. Reeves, J. Yin, Y. Zhu, K. Z. Ibrahim, C. Yang, and V. Vlček, Dynamic mode decomposition for extrapolating nonequilibrium Green’s-function dynamics, Phys. Rev. B 107, 075107 (2023a).
  * Reeves  _et al._ [2023b] C. C. Reeves, Y. Zhu, C. Yang, and V. Vlček, Unimportance of memory for the time nonlocal components of the kadanoff-baym equations, Phys. Rev. B 108, 115152 (2023b).
  * Schlünzen  _et al._ [2017] N. Schlünzen, J.-P. Joost, F. Heidrich-Meisner, and M. Bonitz, Nonequilibrium dynamics in the one-dimensional Fermi-Hubbard model: Comparison of the nonequilibrium Green-functions approach and the density matrix renormalization group method, Phys. Rev. B 95, 165139 (2017).
  * Bonitz  _et al._ [2013] M. Bonitz, S. Hermanns, and K. Balzer, Dynamics of Hubbard Nano-Clusters Following Strong Excitation, Contributions to Plasma Physics 53, 778 (2013).
  * Balzer  _et al._ [2013] K. Balzer, S. Hermanns, and M. Bonitz, The generalized Kadanoff-Baym ansatz. Computing nonlinear response properties of finite systems, Journal of Physics: Conference Series 427, 012006 (2013).
  * Hermanns  _et al._ [2013] S. Hermanns, K. Balzer, and M. Bonitz, Few-particle quantum dynamics–comparing nonequilibrium Green functions with the generalized Kadanoff–Baym ansatz to density operator theory, Journal of Physics: Conference Series 427, 012008 (2013).
  * von Friesen  _et al._ [2010] M. P. von Friesen, C. Verdozzi, and C.-O. Almbladh, Artificial damping in the Kadanoff-Baym dynamics of small Hubbard chains, Journal of Physics: Conference Series 220, 012016 (2010).
  * von Friesen  _et al._ [2009] M. P. von Friesen, C. Verdozzi, and C.-O. Almbladh, Successes and Failures of Kadanoff-Baym Dynamics in Hubbard Nanoclusters, Phys. Rev. Lett. 103, 176404 (2009).
  * [65] See Supplemental Material at  for technical details on each of the methods employed in this paper, as well as the relevant implementation details. Additional results are also included: Comparision of methods in a long range interacting model, asymptotic behaviour of the equilibrium natural occupations with increasing $U$ and analysis of the system entropy with different system sizes.
  * Park and Light [1986] T. J. Park and J. C. Light, Unitary quantum time evolution by iterative Lanczos reduction, J. Chem. Phys. 85, 5870 (1986).
  * Stan  _et al._ [2009] A. Stan, N. Dahlen, and R. van Leeuwen, Time propagation of the Kadanoff–Baym equations for inhomogeneous systems, J. Chem. Phys. 130, 224101 (2009).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Physical Review B 34, 5390 (1986).
  * Joost  _et al._ [2020] J.-P. Joost, N. Schlünzen, and M. Bonitz, G1-G2 scheme: Dramatic acceleration of nonequilibrium Green functions simulations within the Hartree-Fock generalized Kadanoff-Baym ansatz, Phys. Rev. B 101, 245101 (2020).
  * Harsha  _et al._ [2019] G. Harsha, T. M. Henderson, and G. E. Scuseria, Thermofield Theory for Finite-Temperature Coupled Cluster, J. Chem. Theory Comput. 15, 6127 (2019).
  * Yuan  _et al._ [2022] S. Yuan, Y. Chang, and L. K. Wagner, Quantification of electron correlation for approximate quantum calculations, The Journal of Chemical Physics 157, 194101 (2022), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0119260/16552862/194101_1_online.pdf .
  * Reeves and Vlcek [2024] C. Reeves and V. Vlcek, A real-time dyson expansion scheme: Efficient inclusion of dynamical correlations in non-equilibrium spectral propertie (2024), arXiv:2403.07155 [physics.comp-ph] .


Nuclear transparency of the charged hadron produced in the inclusive
electronuclear reaction

Swapan Das 111email<EMAIL_ADDRESS>
Nuclear Physics Division, Bhabha Atomic Research Centre,
Trombay, Mumbai 400085, India
Homi Bhabha National Institute, Anushakti Nagar, Mumbai 400094, India

###### Abstract

The nuclear transparency of the charged hadron produced in the inclusive
$(e,e^{\prime})$ reaction on the nucleus has been calculated using Glauber
model for the nuclear reaction. The color transparency (CT) of the produced
hadron and the short-range correlation (SRC) of the nucleons in the nucleus
have been incorporated in the Glauber model to investigate their effects on
the nuclear transparency of the hadron. The calculated nuclear transparencies
for the proton and pion are compared with the data.

The hadron-nucleus cross section is less than that in the plane wave impulse
approximation (PWIA) because the initial and(or) final state interactions of
the hadron with the nucleus are neglected in PWIA. The difference in the cross
sections can be characterized by the nuclear transparency $T_{A}$, defined [1]
as

$T_{A}=\frac{\sigma_{hA}}{\sigma_{hA(PWIA)}},$ (1)

where $\sigma_{hA}$ represents the hadron-nucleus cross section.

The transverse size $d_{\perp}$ of the hadron produced in the nucleus due to
the space-like high momentum transfer $Q^{2}$ is reduced as $d_{\perp}\sim
1/Q$ [1, 2]. The reduced (in size) hadron is referred as point like
configuration (PLC) [1]. According to Quantum Chromodynamics, a color neutral
PLC has reduced interaction with the nucleon in the nucleus because the sum of
its gluon emission amplitudes cancel [1, 3]. The PLC expands to the size of
physical hadron, as it moves up to a length ($\sim 1$ fm) called hadron
formation length $l_{h}$ [1, 4]:

$l_{h}=\frac{2k_{h}}{\Delta M^{2}},$ (2)

where $k_{h}$ is the momentum of hadron in the laboratory frame. $\Delta
M^{2}$ is related to the mass difference between the hadronic states
originating due to the (anti)quarks fluctuation in PLC. The interaction of PLC
with the nucleon in the nucleus increases, as its size increases during its
passage $l_{h}$ through the nucleus. The decrease in the hadron-nucleon cross
section in the nucleus, as explained by Glauber model [5], leads to the
increase in the hadron-nucleus cross section $\sigma_{hA}$. Therefore, the
transparency $T_{A}$ in Eq. (1) of the hadron raises. The enhancement in
$T_{A}$ due to the above phenomenon is referred as color transparency (CT) of
the hadron. The physics of CT for hadrons have been discussed elaborately in
Refs. [3, 6].

The experiments to search the CT of proton ($p$CT) in the A$(p,pp)$ reactions
done at Brookhaven National Laboratory (BNL) [7] could not confirm the $p$CT.
In fact, the experimental results are not understood [8]. The $p$CT is not
also seen in the A$(e,e^{\prime}p)$ experiment done at Standford Linear
Accelerator Center (SLAC) [9] and Jefferson Laboratory (JLab) [10] for
$0.64\leq Q^{2}\leq 8.1$ GeV2. This experiment done for $8\leq
Q^{2}\mbox{(GeV${}^{2}$)}\leq 14.2$ [11] at the upgraded JLab facility agrees
with the previous observation [9, 10]. Therefore, it appears the PLC required
for CT is unlikely to form for three quarks $(qqq)$ system, such as proton.

Since the meson is a bound state of two quarks (i.e., quark-antiquark) the PLC
formation of it can be more probable than that of the baryon, a three quarks
$(qqq)$ system. The color transparency is unambiguously reported from Fermi
National Accelerator Laboratory (FNAL) [12] in the experiment of the nuclear
diffractive dissociation of pion (of 500 GeV/$c$) to dijets. The color
transparency is also illustrated in the $\pi^{-}$ meson photoproduction [13]
and $\rho^{0}$ meson electroproduction (from nuclei) experiments [14]. Several
authors studied the $\rho$-meson color transparency in the energy region
available at JLab [1, 15].

The nuclear transparency of the $\pi^{+}$ meson produced in the
A$(e,e^{\prime})$ process was measured at JLab for the photon virtuality
$Q^{2}=1.1-4.7$ GeV2 [16]. The data have been understood by the pionic color
transparency ($\pi$CT) [17]. Larson et al., [18] described the momentum
dependence of $\pi$CT in the above reaction. Cosyn et al., [19] studied the
effects of $\pi$CT and nucleon short-range correlation in the pion photo- and
electro- production from nuclei. Larionov et al., [4] estimated the $\pi$CT in
the $(\pi^{-},l^{+}l^{-})$ reaction on nuclei for $p_{\pi}=5-20$ GeV/$c$,
which can be measured at the forth-coming facilities in Japan Proton
Accelerator Research Complex (J-PARC) [20]. This reaction provides the
informations complementary to those obtained from the A$(\gamma^{*},\pi)$
reaction. Miller and Strikman [21] illustrated large CT in the pionic knockout
of the proton off nuclei at the energy 200 GeV available at CERN COMPASS
experiment.

The enhancement in $T_{A}$ due to $\sigma_{hA}$ in Eq. (1) can also occur
because of the short-range correlation (SRC) of nucleon in the nucleus. The
SRC arises because of the repulsive (short-range) interaction between the
nucleons bound in the nucleus. This interaction keeps the bound nucleons apart
($\sim$ 1 fm), which is called nuclear granularity [8]. Therefore, the SRC
prevents the shadowing of the hadron-nucleon interaction due to the
surrounding nucleons present in the nucleus. This occurrence, as elucidated by
Glauber model [5], leads to the enhancement in $\sigma_{hA}$. The SRC is
widely used to investigate various aspects in the nuclear physics [22].

The hadron $h$ is produced in the inclusive A$(e,e^{\prime})$X reaction
because of the interaction of the virtual photon $\gamma^{*}$ (emitted at the
$ee^{\prime}$ vertex) with the nucleus A. In this reaction, the nucleus in the
final state denoted by X is unspecified. The scattering amplitude for the
$\gamma^{*}\mbox{A}\to h\mbox{X}$ transition, according to Glauber model [1],
can be written as

$F_{X0}[({\bf{q-k}}_{h})_{\perp}]=\frac{iq}{2\pi}\int d{\bf
b}e^{i({\bf{q-k}}_{h})_{\perp}\cdot{\bf{b}}}\Gamma^{\gamma^{*}h}_{X0}({\bf
b}),$ (3)

where ${\bf q}$ and ${\bf k}_{h}$ are the momenta of $\gamma^{*}$ and $h$
respectively. $\Gamma^{\gamma^{*}h}_{X0}({\bf b})$ describes the matrix
element for the transition of the nucleus from its initial to final states,
i.e.,

$\Gamma^{\gamma^{*}h}_{X0}({\bf b})=<X|\Gamma^{\gamma^{*}h}_{A}({\bf b},{\bf
r}_{1},...,{\bf r}_{A})|0>,$ (4)

where $|0>$ denotes the ground state of the target nucleus and $|X>$
represents the unspecified nuclear state in the exit channel. The nuclear
profile operator $\Gamma^{\gamma^{*}h}_{A}({\bf b},{\bf r}_{1},...,{\bf
r}_{A})$ [1, 23] is given by

$\Gamma^{\gamma^{*}h}_{A}({\bf b},{\bf r}_{1},...,{\bf
r}_{A})=\sum_{i}\Gamma^{\gamma^{*}h}({\bf b-b}_{i})e^{i({\bf
q-k}_{h})_{\parallel}z_{i}}\Pi^{A-1}_{j\not{=}i}[1-\Gamma^{hN}({\bf
b-b}_{j})\theta(z_{j}-z_{i})].$ (5)

The summation $i$ is taken over the number of nucleons in the nucleus
participated for the hadron production, e.g., the protons in the nucleus take
part to produced charged hadron in the reaction.

$\Gamma^{\gamma^{*}h}({\bf\tilde{b}})$ is the two-body profile function for
the hadron produced from the nucleon, i.e., $\gamma^{*}N\to hN$ process. It is
related to the reaction amplitude $f_{\gamma^{*}h}({\bf\tilde{q}}_{\perp})$
[1] as

$\Gamma^{\gamma^{*}h}({\bf\tilde{b}})=\frac{1}{i2\pi q}\int
d{\bf\tilde{q}}_{\perp}e^{-i{\bf\tilde{q}}_{\perp}\cdot{\bf\tilde{b}}}f_{\gamma^{*}h}({\bf\tilde{q}}_{\perp}).$
(6)

The two-body profile function $\Gamma^{hN}({\bf\tilde{b}})$ is connected to
$hN$ (hadron-nucleon) elastic scattering amplitude
$f_{hN}({\bf\tilde{q}}_{\perp})$ [1, 5] as

$f_{hN}({\bf\tilde{q}^{\prime}}_{\perp})=\frac{ik_{h}}{2\pi}\int
d{\bf\tilde{b}^{\prime}}e^{i{\bf\tilde{q}^{\prime}}_{\perp}\cdot{\bf\tilde{b}^{\prime}}}\Gamma^{hN}({\bf\tilde{b}^{\prime}}).$
(7)

The nuclear states, assuming the independent particle model [24], can be
written in terms of the single particle state $\Phi$ as
$|0>=\Pi^{A}_{l=1}|\Phi_{0}({\bf r}_{l})>$ and $|X>=|\Phi_{X}({\bf
r}_{m})>\Pi^{A-1}_{n{\not=}m}|\Phi_{0}({\bf r}_{n})>$. Using those,
$\Gamma^{\gamma^{*}h}_{X0}({\bf b})$ in Eq. (4) can be written as

$\Gamma^{\gamma^{*}h}_{X0}({\bf b})=\sum_{i}\int d{\bf r}_{i}\Phi^{*}_{X}({\bf
r}_{i})\Gamma^{\gamma^{*}h}({\bf b-b}_{i})e^{i({\bf
q-k}_{h})_{\parallel}z_{i}}\Phi_{0}({\bf r}_{i})D({\bf b},z_{i}),$ (8)

where $D({\bf b},z_{i})$ is given by

$\displaystyle D({\bf b},z_{i})$ $\displaystyle=$
$\displaystyle\Pi^{A-1}_{j{\not=}i}\int d{\bf r}_{j}\Phi^{*}_{0}({\bf
r}_{j})[1-\Gamma^{hN}({\bf b-b}_{j})\theta(z_{j}-z_{i})]\Phi_{0}({\bf r}_{j})$
(9) $\displaystyle=$ $\displaystyle\left[1-\frac{1}{A}\int d{\bf
b}_{j}\Gamma^{hN}({\bf b-b}_{j})\int dz_{j}\theta(z_{j}-z_{i})\varrho({\bf
r}_{j})\right]^{A-1}.$

In this equation, $\varrho({\bf r}_{j})$ is the matter density distribution of
the nucleus, i.e., $\varrho({\bf r}_{j})=A|\Phi_{0}({\bf r}_{j})|^{2}$.
$\varrho({\bf b}_{j},z_{j})$ can be replaced by $\varrho({\bf b},z_{j})$,
since $\Gamma^{hN}({\bf b-b}_{j})$ varies much rapidly than $\varrho({\bf
b}_{j},z_{j})$ [1]. Using Eq. (7) and ${\cal
L}t_{n\to\infty}(1+\frac{x}{n})^{n}=e^{x}$, the above equation can be
simplified to

$D({\bf b},z_{i})\simeq e^{-\frac{1}{2}\sigma^{hN}_{t}[1-i\alpha_{hN}]T({\bf
b},z_{i})},$ (10)

where $\alpha_{hN}$ denotes the ratio of the real to imaginary part of the
hadron-nucleon scattering amplitude $f_{hN}(0)$, and $\sigma^{hN}_{t}$ =
$\frac{4\pi}{k_{h}}Im[f_{hN}(0)]$ is the hadron-nucleon total cross section.
$T({\bf b},z_{i})$ is the partial thickness function of the nucleus, i.e.,

$T({\bf b},z_{i})=\int^{\infty}_{z_{i}}dz_{j}\varrho({\bf b},z_{j}).$ (11)

Using Eq. (8), $F_{X0}[({\bf{q-k}}_{h})_{\perp}]$ in Eq. (3) can be expressed
as

$\displaystyle F_{X0}$ $\displaystyle=$ $\displaystyle\frac{iq}{2\pi}\int
d{\bf b}e^{i({\bf{q-k}}_{h})_{\perp}\cdot{\bf{b}}}\sum_{i}\int d{\bf
r}_{i}\Phi^{*}_{X}({\bf r}_{i})\Gamma^{\gamma^{*}h}({\bf b-b}_{i})e^{i({\bf
q-k}_{h})_{\parallel}z_{i}}\Phi_{0}({\bf r}_{i})D({\bf b},z_{i}),$ (12)
$\displaystyle=$ $\displaystyle\sum_{i}\int d{\bf r}_{i}\Phi^{*}_{X}({\bf
r}_{i})f^{(i)}_{hN}([{\bf{q-k}}_{h}]_{\perp})e^{i({\bf q-k}_{h})\cdot{\bf
r}_{i}}\Phi_{0}({\bf r}_{i})D({\bf r}_{i}).$

$f^{(i)}_{hN}$, defined in Eq. (7), can be considered identically equal for
all nucleons.

The nucleus in the final state $|X>$ differs from its initial state $|0>$
(i.e., ground state) for the charged hadron production, i.e.,
$\Phi_{X}{\not=}\Phi_{0}$ and $F_{00}=0$. To calculate the cross section,
$|F_{X0}|^{2}$ is to multiply by the phase-space of the reaction and that is
to divide by the incident flux. Since the final state $|X>$ of the nucleus is
not detected in the inclusive reaction, the summation over all states $X$ has
to carry out. In the multi-GeV region, the phase space of the reaction can be
considered independent of the state $X$, and therefore, the nuclear
transparency $T_{A}$ can be written [1] as

$T_{A}=\frac{\sum_{X\not{=}0}|F_{X0}|^{2}}{\sum_{X\not{=}0}|F_{X0}|^{2}_{PWIA}}.$
(13)

The hadron-nucleon cross section $\sigma^{hN}_{t}$ in the free-space is used
in Eq. (10) to evaluate $T_{A}$ in Glauber model. To look for the color
transparency (CT), $\sigma^{hN}_{t}$ (according to quantum diffusion model [2,
18]) has to replace by $\sigma^{hN}_{t,CT}$:

$\sigma^{hN}_{t,CT}(Q^{2},l_{z})=\sigma^{hN}_{t}\left[\left\\{\frac{l_{z}}{l_{h}}+\frac{n_{q}^{2}<k^{2}_{t}>}{Q^{2}}\left(1-\frac{l_{z}}{l_{h}}\right)\right\\}\theta(l_{h}-l_{z})+\theta(l_{z}-l_{h})\right],$
(14)

where $Q^{2}$ is the space-like four-momentum transfer, i.e., photon
virtuality. $n_{q}$ denotes the number of valence quak-(anti)quark present in
the hadron, e.g., $n_{q}=2(3)$ for pion (proton) [2]. $k_{t}$ illustrates the
transverse momentum of the (anti)quark: $<k_{t}^{2}>^{1/2}=0.35$ GeV/$c$.
$l_{z}$ is the path length traversed by the hadron after its production. The
hadron formation length $l_{h}(\propto\frac{1}{\Delta M^{2}})$ is already
defined in Eq. (2).

The short-range correlation (SRC) can be incorporated by replacing the nuclear
density distribution $\varrho$ in Eq. (11) by

$\varrho({\bf b},z_{j})\to\varrho({\bf b},z_{j})C(|z_{j}-z_{i}|),$ (15)

where $C(u)$ represents the correlation function [8]. Using the nuclear matter
estimate, it can be written as

$C(u)=\left[1-\frac{h(u)^{2}}{4}\right]^{1/2}[1+f(u)],$ (16)

with $h(u)=3\frac{j_{1}(k_{F}u)}{k_{F}u}$ and $f(u)=-e^{-\alpha u^{2}}(1-\beta
u^{2})$. The Fermi momentum $k_{F}$ is chosen equal to 1.36 fm-1. $C(u)$ with
the parameters $\alpha=1.1$ fm-2 and $\beta=0.68$ fm-2 agrees well that
derived from the many-body calculations [8].

The nuclear transparency $T_{A}$ of the charged hadron, i.e., proton and
$\pi^{+}$ meson, produced in the inclusive electronuclear reaction has been
calculated using Glauber model (GM), where the measured nuclear density
distribution $\varrho(r)$ [25], and hadron-nucleon cross section
$\sigma_{t}^{hN}$ [26] are used. As shown later, the calculated results due to
GM (presented by the dashed curves) underestimate the measured $T_{A}$ for
both proton and pion. Therefore, GM has been modified by taking account of CT
and SRC. Since the CT is energy dependent, the calculated $T_{A}(\pi^{+})$
increases with $Q^{2}$ due to the inclusion of CT in GM. Unlike CT, the SRC is
independent of energy. Therefore, the calculated $T_{A}(\pi^{+})$ due to SRC
added in GM does not illustrate the $Q^{2}$ dependence. The dot-dot-dashed and
dot-dashed curves arise because of the inclusion of CT in GM for $\Delta
M^{2}$, defined in Eq. (2), taken equal to 0.7 and 1.4 GeV2 respectively. The
calculated $T_{A}$ due to SRC incorporated in GM are presented by the solid
curves.

The calculated proton transparency $T_{A}(p)$ vs photon virtuality $Q^{2}$ in
the A$(e,e^{\prime}p)$X reaction is compared with the data in Fig. 1. The data
reported from SLAC [9] and JLab [10, 11] are represented by the white squares
and black circles respectively. Fig. 1(a) shows CT does not exist for the
proton moving through 12C for a wide range of the photon virtuality, i.e.,
$0.64\leq Q^{2}\leq 14.2$ GeV2. This is corroborated by the results for other
nuclei shown in Figs. 1(b) and (c), where the data are available for lesser
range of $Q^{2}$, i.e., $0.64\leq Q^{2}\leq 8.1$ GeV2 for 56Fe and $0.64\leq
Q^{2}\leq 6.77$ GeV2 for 197Au. Therefore, the CT of proton is distinctly
ruled out. Fig. 1 shows the calculated $T_{A}(p)$ due to the inclusion of SRC
in GM reproduce the data reasonably well for all nuclei.

The measured pionic transparency $T_{A}(\pi^{+})$ for $1.1\leq Q^{2}\leq 4.69$
GeV2 in the A$(e,e^{\prime}\pi^{+})$X reaction have been reported from JLab
[16] for 12C, 27Al, 63Cu and 197Au nuclei. The data for all nuclei (except
12C) show the enhancement of $T_{A}(\pi^{+})$ with $Q^{2}$. Proposals are
there to measure $T_{A}(\pi^{+})$ at JLab for higher $Q^{2}$, i.e., $5\leq
Q^{2}\leq 9.5$ GeV2 [3, 27]. Therefore, $T_{A}(\pi^{+})$ for $1.1\leq
Q^{2}\leq 9.5$ GeV2 have been calculated and those are presented in Fig. 2
along with the available data [16]. The calculated results due to GM+CT (i.e.,
$\pi$CT) are accord with both the $Q^{2}$ dependence and magnitude of the
data. This is in concurrence with the earlier calculations [4, 17]. The
calculated $T_{A}(\pi^{+})$ due to GM+SRC do not describe the $Q^{2}$
dependence of the data but the calculated results agree with the large number
of data points within the errors. Therefore, the data of $T_{A}(\pi^{+})$ in
the region of $Q^{2}=5-9.5$ GeV2 are necessary to prove the existence of
$\pi$CT.

The author appreciates Prof. Dipnagkar Dutta for the discussions on the
experimental results, and thanks A. K. Gupta and S. M. Yusuf for their
encouragement to work on theoretical nuclear physics.

## References

  * [1] G. T. Howell and G. A. Miller, Phys. Rev. C 88 (2013) 035202.
  * [2] G. R. Farrar, H. Liu, L. L. Frankfurt and M. I. Strikman, Phys. Rev. Lett. 61 (1988) 686.
  * [3] D. Dutta, K. Hafidi, and M. Strikman, Prog. Part. Nucl. Phys. 69 (2013) 1.
  * [4] A. B. Larionov, M. Strikman and M. Bleicher, Phys. Rev. C 93 (2016) 034618.
  * [5] R. J. Glauber, in Lectures in Theoretical Physics, edited by W. E. Brittin et al. (Interscience, New York, 1959), Vol. I, p. 315; J. M. Eisenberg and D. S. Kolton, Theory of Meson Interaction with Nuclei (John Wiley $\&$ Sons, New York, 1980) p. 158.
  * [6] L. L. Frankfurt, G. A. Miller and M. Strikman, Annu. Rev. Nucl. Part. Sci. 44 (1994) 501; L. Frankfurt and M. Strikman, Phys. Rep. 160 (1988) 235; P. Jain, B. Pire and J. P. Ralston, Phys. Rep. 271 (1996) 67.
  * [7] A. S. Carroll et al., Phys. Rev. Lett. 61 (1988) 1698; I. Mardor et al., Phys. Rev. Lett. 81 (1998) 5085; A. Leksanov et al., Phys. Rev. Lett. 87 (2001) 212301.
  * [8] T.-S. H. Lee and G. A. Miller, Phys. Rev. C 45 (1992) 1863.
  * [9] T. G. O’Neill et al., Phys. Lett. B 351 (1995) 87; N. C. R. Makins et al., Phys. Rev. Lett. 72 (1994) 1986.
  * [10] D. Dutta et al., Phys. Rev. C 68 (2003) 064603; D. Abbott et al., Phys. Rev. Lett 80 (1998) 5072; K. Garrow et al., Phys. Rev. C 66 (2002) 044613.
  * [11] D. Bhetuwal et al.,Phys. Rev. Lett. 126 (2021) 082301.
  * [12] E. M. Aitala et al., Phys. Rev. Lett. 86 (2001) 4773.
  * [13] D. Dutta et al. (E94104 Collaboration), Phys. Rev. C 68 (2003) 021001R.
  * [14] A. Airapetian et al., Phys. Rev. Lett. 90 (2003) 052501; L. El Fassi et al., Phys. Lett. B 712 (2012) 326.
  * [15] B. Z. Kopeliovich, J. Nemchik and I. Schmidt, Phys. Rev. C 76 (2007) 015205; L. Frankfurt, G. A. Miller and M. Strikman, Phys. Rev. C 78 (2008) 015208; K. Gallmeister, M. Kaskulov and U. Mosel, Phys. Rev. C 83 (2011) 015201.
  * [16] B. Clasie et al., Phys. Rev. Lett. 99 (2007) 242502; X. Qian et al., Phys. Rev. C 81 (2010) 055209.
  * [17] M. M. Kaskulov, K. Gallmeister and U. Mosel, Phys. Rev. C 79 (2009) 015207.
  * [18] A. Larson, G. A. Miller and M. Strikman, Phys. Rev. C 74 (2006) 018201.
  * [19] W. Cosyn, M. C. Mart$\i^{\prime}$nez and J. Ryckebusch, Phys. Rev. C 77 (2008) 034602.
  * [20] S. Kumano, in 21st International Symposium on Spin Physics (SPIN 2014) Beijing, China, October 20-24, 2014: arXiv 1504.05264 [hep-ph]; Int. J. Mod. Phys. (Conference Series) 40 (2016) 1660009.
  * [21] G. A. Miller and M. Strikman, Phys. Rev. C 82 (2010) 025205.
  * [22] G. A. Miller and J. E. Spencer, Ann. Phys. (N.Y.) 100 (1976) 562; O. Benhar et al., Phys. Rev. C 44 (1991) 2328; S. Das, Phys. Scr. 96 (2021) 035304.
  * [23] J. H$\ddot{\mbox{u}}$fner, B. Kopeliovich and J. Nemchik, Phys. Lett. B 383 (1996) 362.
  * [24] T. H. Bauer, R. D. Spital, D. R. Yennie and F. M. Pipkin, Rev. Mod. Phys. 50 (1978) 261; Erratum, 51 (1979) 407.
  * [25] C. W. De Jager, H. De Vries and C. De Vries, At. Data and Nucl. Data Tables 14 (1974) 479; 36 (1987) 495.
  * [26] P. A. Zyla et al., (Particle Data Group), Prog. Theor. Exp. 2020 (2020) 083C01; https://pdg.lbl.gov./2020/hadronic-xsections/hadron.html; C. Lechanoine-Leluc and F. Lehar, Rev. Mod. Phys. 65 (1993) 47; D. V. Bugg et al., Phys. Rev. 146 (1966) 980; S. Barshay, C. B. Dover and J. P. Vary, Phys. Rev. C 11 (1975) 360.
  * [27] D. Dutta, private communication.

Figure 1:  (color online). The calculated nuclear transparency of the proton
$T_{A}(p)$ vs. photon virtuality $Q^{2}$. The dashed curve denotes $T_{A}(p)$
calculated using Glauber model (GM). The dot-dot-dashed and dot-dashed curves
illustrate the proton color transparency ($p$CT) for two different values of
$\Delta M^{2}$, see text. The solid curves arise due to the inclusion of the
short-range correlation (SRC) in GM. The data are taken from Refs. [9]-[11].

Figure 2:  (color online). Same as those presented in Fig. 1 but for the pion.
The data are taken from Refs. [16].