# Nonadiabatic Phase Transition with Broken Chiral Symmetry

Bin Yan  Center for Nonlinear Studies, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA  Theoretical Division, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Vladimir Y. Chernyak
Department of Chemistry and Department of Mathematics, Wayne State University,
Detroit, Michigan 48202, USA  Wojciech H. Zurek  Theoretical Division, Los
Alamos National Laboratory, Los Alamos, New Mexico 87545, USA  Nikolai A.
Sinitsyn  Theoretical Division, Los Alamos National Laboratory, Los Alamos,
New Mexico 87545, USA

###### Abstract

We explore nonadiabatic quantum phase transitions in an Ising spin chain with
a linearly time-dependent transverse field and two different spins per unit
cell. Such a spin system passes through critical points with gapless
excitations, which support nonadiabatic transitions. Nevertheless, we find
that the excitations on one of the chain sublattices are suppressed in the
nearly adiabatic regime exponentially. Thus, we reveal a coherent mechanism to
induce exponentially large density separation for different quasiparticles.

The task to separate two slightly different types of particles is often
encountered in both applied and fundamental research. One example is the
apparent asymmetry between matter and antimatter in our Universe. Observations
prove that the symmetry between matter and antimatter was broken, presumably
early on in the history of the Universe. However, it is still unclear how the
subtle $CP$ symmetry violation could lead to the large observed differences at
the cosmological scale Canetti  _et al._ (2012), although it is known that,
when the characteristic asymmetry energy scale is very small, the large
density difference has to be produced during a nonequilibrium process Sakharov
(1991).

The goal of this Letter is to introduce a mechanism for controlling different
quasiparticles separately using exponential sensitivity of quantum
nonadiabatic transitions to symmetry-breaking interactions. Namely, when
parameters of a quantum system change with time adiabatically, the system
remains in the instantaneous eigenstate, e.g., the ground state. However, the
Kibble-Zurek mechanism Kibble (1980, 1976, 2007); Zurek (1985, 1996); del
Campo and Zurek (2014); Keesling  _et al._ (2019); Rams  _et al._ (2019); Qiu
_et al._ (2020) predicts that the number of nonadiabatic excitations is not
suppressed exponentially if a macroscopic system passes through a quantum
critical point at which the energy gap to excitations closes. Without an
asymmetry of interactions, different particle species would pass through the
same phase transition simultaneously. However, even a small asymmetry
generally opens the gap to certain excitation types, even though the critical
point is not destroyed. Thus, we can harvest the excitations of one type and
suppress the others.

A broadly known quantum example that confirms the Kibble-Zurek mechanism is
the model of a uniform Ising spin chain in a transverse magnetic field
Dziarmaga (2005); Zurek  _et al._ (2005); Heyl  _et al._ (2013), with the
Hamiltonian

$\hat{H}_{u}=\sum_{n=1}^{N}\left[B\hat{\sigma}_{n}^{z}+J\hat{\sigma}_{n}^{x}\hat{\sigma}_{n+1}^{x}\right],\quad
B=-\beta t,$ (1)

where $\hat{\sigma}_{n}^{x,z}$ are Pauli operators for the $n$th spin, $J$ is
the spin-spin coupling, and $B$ is the transverse time-dependent field that
changes with rate $\beta>0$. This model has two Dirac points at $B=\pm J$, at
which its spectral gaps close and which mark boundaries between three phases.
The phase with strong quantum correlations (Fig. 1, top) contains the point
$B=0$ with two degenerate ground states:

$|\rightarrow\,\rightarrow\,\cdots\rightarrow\,\rightarrow\rangle\quad{\rm
and}\quad|\leftarrow\,\leftarrow\,\cdots\leftarrow\,\leftarrow\rangle.$ (2)

Figure 1: Top: phase diagram of Ising chain (1) in a transverse magnetic
field. Dashed horizontal lines mark the critical points with the spectral gap
closings. At $B=\pm\infty$, the ground state is fully spin polarized along the
$z$ axis, whereas at $B=0$ the ground state is a superposition of
$x$-polarized spin states. Bottom: spin chain model (5) with broken chiral
symmetry. We consider the case of $b_{1}>b_{2}$. In the adiabatic limit, the
initial spin polarized ground state (top) at $t=-\infty$ transfers to the spin
polarized state at $t=+\infty$ with all spins flipped. The nonadiabatic
excitations (nonflipped spins) on the odd-site sublattice (red) are suppressed
as a power law (31), respecting the Kibble-Zurek mechanism. In contrast,
nonadiabatic excitations on the even-site sublattice (blue) are suppressed
exponentially (25).

For $B(t)=-\beta t$, the model is exactly solvable Dziarmaga (2005). If the
system starts from the fully polarized ground state

$|\downarrow\,\downarrow\,\downarrow\,\downarrow\ldots\downarrow\,\downarrow\rangle$
(3)

as $t\rightarrow-\infty$, then in the thermodynamic limit of this solution,
after the passage through a critical point, the number of nonadiabatic
excitations follows a power law:

$\rho\sim J^{-1}\beta^{1/2}.$ (4)

Let us now explore robustness of the Ising model predictions by considering a
spin chain whose unit cell contains two different spins, with the Hamiltonian

$\displaystyle\hat{H}=$
$\displaystyle\sum_{n=1}^{N}\left[-\frac{b_{1}t}{2}\hat{\sigma}_{2n}^{z}-\frac{b_{2}t}{2}\hat{\sigma}_{2n+1}^{z}\right]$
(5)
$\displaystyle+\sum_{n=1}^{N}\left[J_{1}\hat{\sigma}_{2n}^{x}\hat{\sigma}_{2n+1}^{x}+J_{2}\hat{\sigma}_{2n+1}^{x}\hat{\sigma}_{2n+2}^{x}\right].$

The difference between $b_{1}$ and $b_{2}$ is due to different spin $g$
factors, and $J_{1}\neq J_{2}$ reflects the lack of inversion symmetry at zero
field. A periodic boundary condition is imposed. Figure 1 (bottom) illustrates
the structure of this spin chain.

Without loss of generality, here and in what follows, we assume $b_{1}>b_{2}$;
namely, the spins on the even-site sublattice have a larger $g$ factor than
those on the odd sites. We note that, for weak symmetry breaking,
$|J_{1}-J_{2}|\ll|J_{1}+J_{2}|$ and $|b_{1}-b_{2}|\ll|b_{1}+b_{2}|$, the
ground state polarizations at $B=0,\pm\infty$ are the same for spins with odd
and even indices, and the spin excitations on different sublattices then
resemble particles of two slightly different types of matter.

Via the Jordan-Wigner transformation, Hamiltonian (5) reduces to a quadratic
form of fermionic operators, $\hat{c}$ and $\hat{d}$, on the even and odd
sites, respectively. The details are presented in Appendix A of Supplemental
Material (SM). Thus,

$\hat{\sigma}_{2n}^{z}=1-2\hat{c}^{\dagger}_{n}\hat{c}_{n},\ \ \
\hat{\sigma}_{2n+1}^{z}=1-2\hat{d}^{\dagger}_{n}\hat{d}_{n}.$ (6)

It is convenient to work with the Fourier transformed operators

$\hat{c}_{p}=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\hat{c}_{n}e^{-ipn},\ \ \
\hat{c}^{\dagger}_{p}=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\hat{c}^{\dagger}_{n}e^{ipn},$
(7)

and similarly defined $\hat{d}_{p}$ and $\hat{d}^{\dagger}_{p}$. We will
assume that $N$ is even, so the momentum takes discrete values
$p=\pm(2k-1)\pi/N$, $k=1,\ldots,N/2$. Hamiltonian (5) then is

$\hat{H}=\sum_{p>0}\hat{\bm{a}}^{\dagger}_{p}{H}_{p}\hat{\bm{a}}_{p},$ (8)

where $\hat{\bm{a}}_{p}$ and ${H}_{p}$ are given by, respectively,

$\hat{\bm{a}}_{p}\equiv\begin{pmatrix}\hat{c}_{p}\\\ \hat{c}^{\dagger}_{-p}\\\
\hat{d}_{p}\\\ \hat{d}^{\dagger}_{-p}\\\ \end{pmatrix},\ \ \
H_{p}\equiv\begin{pmatrix}b_{1}t&0&g&\gamma\\\ 0&-b_{1}t&-\gamma&-g\\\
g^{*}&-\gamma^{*}&b_{2}t&0\\\ \gamma^{*}&-g^{*}&0&-b_{2}t\\\ \end{pmatrix},$
(9)

with the couplings

$g\equiv J_{1}+J_{2}e^{-ip},\ \ \gamma\equiv J_{1}-J_{2}e^{-ip}.$ (10)

In the Heisenberg picture, the evolution of $\hat{\bm{a}}_{p}$ is described by
a Schrödinger-like equation:

$i\hat{\bm{a}}_{p}(t)=H_{p}(t)\hat{\bm{a}}_{p}(t),$ (11)

and the initial ground state (3) corresponds to the initially fully filled
Fermi sea, for all fermions as $t\rightarrow-\infty$.

Figure 2: The time-dependent spectrum of the Hamiltonian $H_{p}(t)$ from Eq.
(9). (a) Two exact level crossing points appear at $p=0$. (b) At finite $p$,
which is here $p=0.3$, small gaps open up near places of the former crossing
points. Other parameters are $b_{1}=2$, $b_{2}=1$, $J_{1}=3$, and $J_{2}=1$ in
both figures. The numbers $1$-$4$ in (a) mark the diagonal elements of
$H_{p}$, with which the corresponding eigenenergies merge as
$t\rightarrow\pm\infty$, e.g., $1$ is for $(H_{p})_{11}$. The dashed blue and
green arrows in (b) show two semiclassical trajectories that contribute to the
amplitude $S_{33}$. Each trajectory jumps through one of the gaps but avoids
the other gap.

In Fig. 2, we show the time-dependent eigenvalues of the matrix $H_{p}$ for a
broken chiral symmetry: $J_{1}\neq J_{2}$ and $b_{1}>b_{2}$. Figure 2(a) shows
that at $p=0$ two of the four energy levels experience exact level crossings,
whereas Fig. 2(b) illustrates that this degeneracy is lifted for nonzero $p$.
This means that the chiral asymmetry does not destroy the Dirac points in the
spectrum, which must produce a power-law density of excitations according to
the Kibble-Zurek mechanism.

We will consider the evolution with the Hamiltonian (5) during the time
interval $t\in(-\infty,+\infty)$. In terms of the phase diagram in Fig. 1,
this means that the system passes through two phase transition points. In the
adiabatic limit, the initial ground state (3) transfers, as
$t\rightarrow+\infty$, into the ground state with all spins up

$|\uparrow\,\uparrow\,\uparrow\,\uparrow\ldots\uparrow\,\uparrow\rangle.$ (12)

This is the fermionic vacuum, and the free fermions are the elementary
excitations. The numbers of excitations (nonflipped spins) on the even and odd
sublattices at the end of the protocol are given by, respectively,

$\displaystyle N_{\rm
ex}^{e}\equiv\sum_{n=1}^{N}\langle\hat{c}_{n}^{\dagger}\hat{c}_{n}\rangle=\sum_{p}^{N}\langle\hat{c}_{p}^{\dagger}\hat{c}_{p}\rangle,$
(13) $\displaystyle N_{\rm
ex}^{o}\equiv\sum_{n=1}^{N}\langle\hat{d}_{n}^{\dagger}\hat{d}_{n}\rangle=\sum_{p}^{N}\langle\hat{d}_{p}^{\dagger}\hat{d}_{p}\rangle,$
(14)

where, in the Schrödinger picture, the averaging is taken over the state at
$t=+\infty$.

The Hamiltonian (9) has the form of a multistate Landau-Zener (MLZ) model,
i.e., can be written as $H(t)=A+Bt$ with Hermitian matrices $A$ and $B$. The
MLZ models have been extensively studied previously Brundobler and Elser
(1993); Sinitsyn  _et al._ (2018); Sinitsyn (2015). Our case with $J_{1}\neq
J_{2}$ is generally not solvable Malla and Raikh (2017) but the MLZ theory
provides exact formulas for some of the evolution matrix elements Li  _et al._
(2017), which we will utilize.

Let us define the $S$ matrix

$S=S(p)\equiv U^{p}(T,-T)_{T\rightarrow\infty},$ (15)

where $U^{p}(T,-T)$ is the evolution matrix over the time interval
$t\in(-T,T)$ with the Hamiltonian ${H}_{p}$. We can say that the operators at
$t=\pm\infty$ are related by

$\hat{\bm{a}}_{p}(+\infty)=S\hat{\bm{a}}_{p},$ (16)

where $\hat{\bm{a}}_{p}\equiv\hat{\bm{a}}_{p}(-\infty)$. Hence,

$\displaystyle\hat{c}_{p}(+\infty)=S_{11}\hat{c}_{p}+S_{12}\hat{c}^{\dagger}_{-p}+S_{13}\hat{d}_{p}+S_{14}\hat{d}^{\dagger}_{-p},$
(17)
$\displaystyle\hat{c}^{\dagger}_{p}(+\infty)=S^{*}_{11}\hat{c}^{\dagger}_{p}+S^{*}_{12}\hat{c}_{-p}+S^{*}_{13}\hat{d}^{\dagger}_{p}+S^{*}_{14}\hat{d}_{-p}.$

The number of excitations (13) can be evaluated in the Heisenberg picture,
i.e., taking the average of the operators at $t=+\infty$ with respect to the
initial state (3):

$N_{\rm ex}^{e}(+\infty)=\sum_{p}\ (|S_{11}|^{2}+|S_{13}|^{2}).$ (18)

The survival amplitudes for states with the highest-energy level slopes are
known exactly for any MLZ model Brundobler and Elser (1993):

$S_{11}=S_{22}=e^{-\pi|g|^{2}/(b_{1}-b_{2})-\pi|\gamma|^{2}/(b_{1}+b_{2})}.$
(19)

Another exact relation of MLZ theory is for the transition amplitudes between
the levels with the two highest slopes Li  _et al._ (2017):

$S_{11}S_{33}+|S_{13}|^{2}=e^{-2\pi|\gamma|^{2}/(b_{1}+b_{2})}.$ (20)

This does not fix $S_{13}$ because $S_{33}$ is not known. Fortunately, model
(9) has discrete symmetries leading to

$\displaystyle S_{11}+S_{22}=S_{33}+S_{44},$ (21) $\displaystyle
S_{33}(-p)=S_{44}(p),$

the derivations of which are given in Appendix B of SM. Equations (19)-(21)
then predict

$\sum_{p}|S_{13}|^{2}=\sum_{p}\left[e^{-2\pi|\gamma|^{2}/(b_{1}+b_{2})}-|S_{11}|^{2}\right].$
(22)

This formula is exact. Hence, without approximations

$N^{e}_{\rm ex}(+\infty)=\sum_{p}\ e^{-2\pi|\gamma|^{2}/(b_{1}+b_{2})}.$ (23)

Figure 3: (a) Density of excitations of the even-site sublattice as a function
of the quenching rate, $b_{1}=2\beta$, $b_{2}=\beta$. The black solid curve
and blue dashed curve correspond to, respectively, the exact prediction (24)
at $J_{1}=J_{2}=1.5$ and perturbed coupling $J_{1}=2$, $J_{2}=1$. (b) The
numerically calculated transition probability $\sum_{p}|S_{33}|^{2}$ at
$J_{1}=2$, $J_{2}=1$, $b_{1}=2\beta$, and $b_{2}=\beta$. Triangles are the
numerical data. The black dashed line is the power-law scaling (31)
prediction. The insets resolve the fast Stueckelberg oscillations around the
mean power-law tail. The filled triangles are the same points in the main
figure and the insets.

As $N\rightarrow\infty$, we replace
$\sum_{p}/N\rightarrow\int_{-\pi}^{\pi}dp/(2\pi)$ and find the density of
excitations (per unit cell) in the thermodynamic limit

$\rho^{e}_{\rm
ex}(+\infty)=\exp\left(-\frac{2\pi(J_{1}^{2}+J_{2}^{2})}{|b_{1}+b_{2}|}\right){I}_{0}\left(\frac{4\pi
J_{1}J_{2}}{b_{1}+b_{2}}\right),$ (24)

where ${I}_{0}(x)$ is the modified Bessel function of the first kind. In the
leading order of small $b_{1}+b_{2}$

$\rho^{e}_{\rm
ex}(+\infty)\approx\exp{\left(-\frac{2\pi(J_{1}-J_{2})^{2}}{|b_{1}+b_{2}|}\right)}\frac{\sqrt{b_{1}+b_{2}}}{2\pi\sqrt{2J_{1}J_{2}}}.$
(25)

This is the central result of our Letter. For any $J_{1}\neq J_{2}$, the
excitations on the even sublattice are suppressed exponentially rather than by
a power law. The latter is found only for a symmetric case with $J_{1}=J_{2}$
but $b_{1}>b_{1}$. Complete exact solution for this special case is presented
in Appendix D of SM. We also note that the above solutions do not smoothly
carry over the $b_{1}=b_{2}$ point, to which our exact MLZ analysis does not
apply because of the degeneracy of the time-dependent Zeeman coupling.

Equation (25) does not contradict the Kibble-Zurek mechanism, because the
power law for excitation density does appear on the odd-site sublattice.
Consider

$N^{o}_{\rm ex}(+\infty)=\sum_{p}(|S_{33}|^{2}+|S_{31}|^{2}).$ (26)

Here, $|S_{31}|^{2}$ can be obtained exactly, because $|S_{31}|=|S_{13}|$.
According to Eq. (22), it produces an exponentially suppressed contribution.
We also prove in Appendix C of SM that
$\sum_{p}|S_{33}|^{2}\geq\sum_{p}|S_{11}|^{2}$, so $N^{o}_{\rm ex}\geq
N^{e}_{\rm ex}$. Hence, only $S_{33}$ contains the information about the
power-law scaling.

The exact crossings of the second and third adiabatic levels of $H_{p}$ happen
at $p=0$ and time moments

$t_{1,2}=\mp 2\sqrt{J_{1}J_{2}/(b_{1}b_{2})}.$ (27)

We interpret $t_{1,2}$ as the moments of the phase transitions. By setting
$t=t_{1,2}+\delta t$ and $p=\delta p$, with small $\delta p$ and $\delta t$,
we find an effective Hamiltonian for the interactions within the subspace of
$H_{p}$ spanned by the two eigenstates, whose energies become nearly
degenerate near the critical points:

$\displaystyle\hat{H}_{\rm
eff}\approx\frac{2\sqrt{J_{1}J_{2}b_{1}b_{2}}(\sqrt{J_{1}J_{2}}\delta
p\hat{\tau}_{x}+\sqrt{b_{1}b_{2}}\delta
t\hat{\tau_{z}})}{\sqrt{(b_{2}J_{1}+b_{1}J_{2})(b_{1}J_{1}+b_{2}J_{2})}},$
(28)

where $\hat{\tau}_{x}$ and $\hat{\tau}_{z}$ are the Pauli operators that act
in subspace of two states with closest to one Dirac point energies. The
Landau-Zener formula applied to Hamiltonian (28) returns the probability of a
nonadiabatic excitation after a transition through one Dirac point:

$P_{\rm ex}=\exp\left({-\frac{2\pi(\delta
p)^{2}(J_{1}J_{2})^{3/2}}{\sqrt{(b_{2}J_{1}+b_{1}J_{2})(b_{1}J_{1}+b_{2}J_{2})}}}\right).$
(29)

Our system passes through two Dirac points. An excitation created near the
first point should remain excited after the second point. The system may also
not produce an excitation during the first crossing but produce it at the
second one. Generally, there is interference between such evolution paths
(Stueckelberg’s oscillations). If we disregard this interference, the
probabilities of the two possibilities simply sum, i.e.,

$P_{\rm 2ex}=2P_{\rm ex}(1-P_{\rm ex}).$ (30)

Integrating $P_{\rm 2ex}$ over $\delta p\in(-\infty,+\infty)$, we find that in
the adiabatic limit

$\sum_{p}\frac{|S_{33}|^{2}}{N}\approx\frac{2-\sqrt{2}}{4\pi(J_{1}J_{2})^{3/4}}\left[(b_{2}J_{1}+b_{1}J_{2})(b_{1}J_{1}+b_{2}J_{2})\right]^{1/4}.$
(31)

This is the estimate of the excitation density on the odd-site sublattice in
the nearly adiabatic limit. In Fig. 3, the numerical simulations confirm the
power-law scaling trend for the parameter dependence in Eq. (31), which is
modulated by fast Stueckelberg oscillations. We also note that Eq. (31)
reproduces correctly the scaling prediction of the uniform chain model (4) if
we set $b_{1}=b_{2}=2\beta$ and $J_{1}=J_{2}$.

Hence, as $t\rightarrow+\infty$, the remaining excitations on the odd-site
sublattice are suppressed according to a power law, in contrast to the
exponential suppression on the even-site sublattice. At the intermediate
moment $t=0$, excitations take the form of superpositions of kinks. The
asymmetry appears then, too, i.e., some of the excitations are exponentially
suppressed in comparison to the others but this effect does not have a simple
interpretation in terms of different kink types. We leave such details to
Appendix E of SM.

Our results illustrate how the Kibble-Zurek mechanism works when a part of a
system is not observable. The latter may happen in atomic gases, in which
optically observable spins interact via spins of other atomic species Roy  _et
al._ (2015). In the asymmetric Ising chain, if the spin excitations on the odd
sites are not observable, the even-site spins still create a visible
ferromagnetic state at zero external field and polarize when the external
field is strong. Hence, looking only at the even spins, one can expect that
this system goes through the same phase transitions as the uniform Ising chain
and, thus, produces excitations with a power-law scaling. Our solution shows,
however, that all observable spin excitations are suppressed in this case
exponentially. A power-law tail of the excitations is hidden then in the
nonobservable spins.

In summary, we demonstrated that small differences between interacting spins
in a simple spin chain can lead to an exponentially large effect when passing
slowly through a phase transition. The underlying mechanism arguably processes
universalities: Strong but symmetric interactions compete with each other near
a critical point, where subtle interaction differences play a decisive role.
As long as a symmetry breaking opens the gap for certain types of
quasiparticles while not destroying the critical point, the exponentially
large density separation of different quasiparticles should happen after the
passage through the phase transition independently of the model’s microscopic
details.

The Kibble-Zurek mechanism for Ising-type quantum phase transition has been
recently studied experimentally in ultracold atoms Keesling  _et al._ (2019).
In such systems, an asymmetry can be introduced by placing atoms in a periodic
potential without the inversion symmetry. If the time-linear field ramp is
induced by changing the ac frequency in the rotating-wave approximation, the
sweeping of the detuning frequency across an optical resonance effectively
mimics the field changes in range $B_{z}\in(+\infty,-\infty)$ Sun and Sinitsyn
(2015); Brierley  _et al._ (2012). Hence, the demonstration of our effect
requires only simple modification of the control protocols in already
accessible for studies dynamic phase transitions.

Let us also speculate about possible applications. Different isotopes have
different spin interactions in an ultracold atomic mixture Tomza (2015).
Hence, an adiabatic passage through a quantum critical point can induce spin
excitations in this mixture, overwhelmingly, for only one of the isotopes. One
can then separate such excited atoms using the magnetic deflection approach
from Ref. Mazur  _et al._ (2014), and, thus, develop a technology for isotope
separation. The asymmetry of matter and antimatter in our Universe may also be
viewed now as a result of a hypothetical transition through a quantum critical
point during cosmic inflation, when the matter could be excited from the
vacuum via quantum nonadiabatic processes Fedderke  _et al._ (2015); Enomoto
and Matsuda .

The dynamic phase transitions are found broadly, from cosmological models to
the experiments with superconductors and ultracold atoms. They are likely
responsible for the scaling of mistakes that limit the quantum annealing
computation techniques Pearson  _et al._ (2019). Fortunately, different
systems show universal behavior, which is driven by relatively simple effects
near the critical point. The quasiparticle separation is one of such effects
that remain relevant near the critical point and, thus, can be used for
control of complex quantum systems.

###### Acknowledgements.

This work was carried out under the support of the U.S. Department of Energy,
Office of Science, Basic Energy Sciences, Materials Sciences and Engineering
Division, Condensed Matter Theory Program. B.Y. also acknowledges partial
support from the Center for Nonlinear Studies.

## References

  * Canetti  _et al._ (2012) L. Canetti, M. Drewes, and M. Shaposhnikov, “Matter and antimatter in the universe,” New Journal of Physics 14, 095012 (2012).
  * Sakharov (1991) A. D. Sakharov, “Violation of CP invariance, C asymmetry, and baryon asymmetry of the universe,” Soviet Physics Uspekhi 34, 392 (1991).
  * Kibble (1980) T. W. B. Kibble, “Some implications of a cosmological phase transition,” Physics Reports 67, 183–199 (1980).
  * Kibble (1976) T. W. B. Kibble, “Topology of cosmic domains and strings,” Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Kibble (2007) T. W. B. Kibble, “Phase-transition dynamics in the lab and the universe,” Physics Today 60, 47–52 (2007).
  * Zurek (1985) W. H. Zurek, “Cosmological experiments in superfluid helium?” Nature 317, 505–508 (1985).
  * Zurek (1996) W. H. Zurek, “Cosmological experiments in condensed matter systems,” Physics Reports 276, 177–221 (1996).
  * del Campo and Zurek (2014) A. del Campo and W. H. Zurek, “Universality of phase transition dynamics: Topological defects from symmetry breaking,” International Journal of Modern Physics. A, Particles and Fields, Gravitation, Cosmology 29, 1430018 (2014).
  * Keesling  _et al._ (2019) A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, “Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator,” Nature 568, 207–211 (2019).
  * Rams  _et al._ (2019) M. M. Rams, J. Dziarmaga, and W. H. Zurek, “Symmetry breaking bias and the dynamics of a quantum phase transition,” Physical Review Letters 123, 130603 (2019).
  * Qiu  _et al._ (2020) L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, H.-X. Yang, T. Tian, Y. Xu, and L.-M. Duan, “Observation of generalized Kibble-Zurek mechanism across a first-order quantum phase transition in a spinor condensate,” Science Advances 6, eaba7292 (2020).
  * Dziarmaga (2005) J. Dziarmaga, “Dynamics of a quantum phase transition: exact solution of the quantum ising model,” Physical Review Letters 95, 245701 (2005).
  * Zurek  _et al._ (2005) W. H. Zurek, U. Dorner, and P. Zoller, “Dynamics of a quantum phase transition,” Physical Review Letters 95, 105701 (2005).
  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, “Dynamical quantum phase transitions in the transverse-field ising model,” Physical Review Letters 110, 135704 (2013).
  * Brundobler and Elser (1993) S. Brundobler and V. Elser, “S-matrix for generalized Landau-Zener problem,” Journal of Physics A: Mathematical and General 26, 1211 (1993).
  * Sinitsyn  _et al._ (2018) N. A. Sinitsyn, E. A. Yuzbashyan, V. Y. Chernyak, A. Patra, and C. Sun, “Integrable Time-Dependent quantum hamiltonians,” Physical Review Letters 120, 190402 (2018).
  * Sinitsyn (2015) N. A. Sinitsyn, “Solvable four-state Landau-Zener model of two interacting qubits with path interference,” Physical Review B 92, 205431 (2015).
  * Malla and Raikh (2017) R. K. Malla and M. E. Raikh, “Loss of adiabaticity with increasing tunneling gap in nonintegrable multistate Landau-Zener models,” Physical Review B 96, 115437 (2017).
  * Li  _et al._ (2017) F. Li, C. Sun, V. Y. Chernyak, and N. A Sinitsyn, “Multistate Landau-Zener models with all levels crossing at one point,” Physical Review A 96, 022107 (2017).
  * Roy  _et al._ (2015) D. Roy, L. Yang, S. A. Crooker, and N. A. Sinitsyn, “Cross-correlation spin noise spectroscopy of heterogeneous interacting spin systems,” Scientific Reports 5, 9573 (2015).
  * Sun and Sinitsyn (2015) C. Sun and N. A. Sinitsyn, “Exact transition probabilities for a linear sweep through a Kramers-Kronig resonance,” Journal of Physics A: Mathematical and Theoretical 48, 505202 (2015).
  * Brierley  _et al._ (2012) R. T. Brierley, C. Creatore, P. B. Littlewood, and P. R. Eastham, “Adiabatic state preparation of interacting two-level systems,” Physical Review Letter 109, 043002 (2012).
  * Tomza (2015) M. Tomza, “Energetics and control of ultracold Isotope-Exchange reactions between heteronuclear dimers in external fields,” Physical Review Letters 115, 063201 (2015).
  * Mazur  _et al._ (2014) T. R. Mazur, B. Klappauf, and M. G. Raizen, “Demonstration of magnetically activated and guided isotope separation,” Nature Physics 10, 601–605 (2014).
  * Fedderke  _et al._ (2015) M. A. Fedderke, E. W. Kolb, and M. Wyman, “Irruption of massive particle species during inflation,” Physical Review D 91, 063505 (2015).
  * (26) S. Enomoto and T. Matsuda, “The exact WKB for cosmological particle production,”  arXiv:2010.14835 (2020) .
  * Pearson  _et al._ (2019) A. Pearson, A. Mishra, Hen I., and D. A. Lidar, “Analog errors in quantum annealing: doom and hope,” NPJ Quantum Information 5 (2019).

Supplemental Material for
“Nonadiabatic Phase Transition with Broken Chiral Symmetry”

##  Appendix A Jordan-Wigner transformation

For the model Hamiltonian

$\hat{H}=\sum_{n=1}^{N}\left[-\frac{b_{1}t}{2}\hat{\sigma}_{2n}^{z}-\frac{b_{2}t}{2}\hat{\sigma}_{2n+1}^{z}\right]+\sum_{n=1}^{N}\left[J_{1}\hat{\sigma}_{2n}^{x}\hat{\sigma}_{2n+1}^{x}+J_{2}\hat{\sigma}_{2n+1}^{x}\hat{\sigma}_{2n+2}^{x}\right]$
(32)

with periodic boundary condition, we define two seperate sets of Jordan-Wigner
transformation for the sublattices with even and odd indices, i.e.,

$\displaystyle\hat{\sigma}_{2n}^{z}=1-2\hat{c}^{\dagger}_{n}\hat{c}_{n},$ (33)
$\displaystyle\hat{\sigma}_{2n+1}^{z}=1-2\hat{d}^{\dagger}_{n}\hat{d}_{n},$
$\displaystyle\hat{\sigma}_{2n}^{x}=-(\hat{c}_{n}^{\dagger}+\hat{c}_{n})\prod_{m<n}(1-2\hat{c}^{\dagger}_{m}\hat{c}_{m})(1-2\hat{d}^{\dagger}_{m}\hat{d}_{m}),$
$\displaystyle\hat{\sigma}_{2n+1}^{x}=-(\hat{d}_{n}^{\dagger}+\hat{d}_{n})(1-2\hat{c}_{n}^{\dagger}\hat{c}_{n})\prod_{m<n}(1-2\hat{c}^{\dagger}_{m}\hat{c}_{m})(1-2\hat{d}^{\dagger}_{m}\hat{d}_{m}),$

where

$\displaystyle\\{\hat{c}_{n},\hat{c}_{n^{\prime}}\\}=\\{\hat{c}^{\dagger}_{n},\hat{c}^{\dagger}_{n^{\prime}}\\}=\\{\hat{d}_{n},\hat{d}_{n^{\prime}}\\}=\\{\hat{d}^{\dagger}_{n},\hat{d}^{\dagger}_{n^{\prime}}\\}=0,$
$\displaystyle\\{\hat{c}_{n},\hat{d}_{n^{\prime}}\\}=\\{\hat{c}^{\dagger}_{n},\hat{d}^{\dagger}_{n^{\prime}}\\}=\\{\hat{c}_{n},\hat{d}^{\dagger}_{n^{\prime}}\\}=\\{\hat{c}^{\dagger}_{n},\hat{d}_{n^{\prime}}\\}=0,$
$\displaystyle\\{\hat{c}_{n},\hat{c}^{\dagger}_{n^{\prime}}\\}=\\{\hat{d}_{n},\hat{d}^{\dagger}_{n^{\prime}}\\}=\delta_{n,n^{\prime}}.$

Their Fourier transforms are defined as

$\hat{c}_{p}=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\hat{c}_{n}e^{-ipn},\ \ \
\hat{c}^{\dagger}_{p}=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\hat{c}^{\dagger}_{n}e^{ipn}.$
(34)

Due to the periodic boundary condition of the spin chain, the discrete
momentum takes values in the range $(-\pi,\pi)$,

$\displaystyle p=\pm(2k-1)\pi/N,\ \ k=1,\ldots,N/2,\ \ N\text{ even};$ (35)
$\displaystyle p=\pm 2k\pi/N,\ \ k=1,\ldots,(N-1)/2,\ \ N\text{ odd}.$

In terms of the fermionic operators, the original spin chain Hamiltonian (32)
transforms to

$\hat{H}=\sum_{p>0}\hat{\bm{a}}^{\dagger}_{p}H_{p}\hat{\bm{a}}_{p},$ (36)

where $\hat{\bm{a}}_{p}$ and $H_{p}$ are defined as

$\hat{\bm{a}}_{p}\equiv\begin{pmatrix}\hat{c}_{p}\\\ \hat{c}^{\dagger}_{-p}\\\
\hat{d}_{p}\\\ \hat{d}^{\dagger}_{-p}\\\ \end{pmatrix},\ \ \
H_{p}\equiv\begin{pmatrix}b_{1}t&0&g&\gamma\\\ 0&-b_{1}t&-\gamma&-g\\\
g^{*}&-\gamma^{*}&b_{2}t&0\\\ \gamma^{*}&-g^{*}&0&-b_{2}t\\\ \end{pmatrix},$
(37)

with the couplings

$g\equiv J_{1}+J_{2}e^{-ip},\ \ \gamma\equiv J_{1}-J_{2}e^{-ip}.$ (38)

##  Appendix B Exact solution for $\sum_{p}|S_{13}|^{2}$

Consider the $S$ matrix in the 4-level Landau-Zener problem of the Hamiltonian
$H_{p}$, which is defined in Eq. (37). Here, the $S$ matrix refers to the
unitary matrix

$S=S(p)\equiv U^{p}(T,-T)_{T\rightarrow\infty},$ (39)

where $U^{p}(T,-T)$ is the evolution matrix over the time interval
$t\in(-T,T)$ with the Hamiltonian $H_{p}$.

This model belongs to a particularly well-understood subclass of multilevel
Landau-Zener (MLZ) problems. Namely, it contains bipartite interactions and
all diabatic levels are crossing in one point in the time-energy diagram Li
_et al._ (2017). Here “bipartite” means that all states can be partitioned
into two groups, such that states in one group are coupled directly only to
states in the other group; and the diabatic levels are linearly time-dependent
diagonal elements of the Hamiltonian. Following Ref. Li  _et al._ (2017), we
list here three properties of the $S$ matrix that are essential for solving
the exact solutions:

(i) As in all other MLZ systems, the amplitudes of surviving on the levels
with the lowest and highest slopes are given by the Brundobler-Elser formula:

$S_{11}=S_{22}=e^{-\frac{\pi|g|^{2}}{b_{1}-b_{2}}-\frac{\pi|\gamma|^{2}}{b_{1}+b_{2}}}.$
(40a) Considering the definition of $g$ and $\gamma$ in (38), this means also
that $S_{11}(p)=S_{22}(p)=S_{11}(-p)=S_{22}(-p).$ (40b)

(ii) The other two diagonal elements of the scattering matrix, $S_{33}$ and
$S_{44}$, can be different from each other and from Eq. (40a), but they are
real and satisfying the relations (see Eq. (21) in Ref. Li  _et al._ (2017)):

$S_{11}S_{33}+|S_{13}|^{2}=e^{-\frac{2\pi|\gamma|^{2}}{b_{1}+b_{2}}}$ (41a)
$S_{44}S_{22}+|S_{42}|^{2}=e^{-\frac{2\pi|\gamma|^{2}}{b_{1}+b_{2}}}$ (41b)

(iii) Following Ref. Li  _et al._ (2017), we can derive one more constraint
for the diagonal elements of the $S$ matrix. Define a matrix

$\Theta\equiv\begin{pmatrix}\hat{\mathbb{I}}_{2}&0\\\
0&-\hat{\mathbb{I}}_{2}\\\ \end{pmatrix},$

where $\hat{\mathbb{I}}_{2}$ is the $2\times 2$ identity matrix. It can be
verified that Hamiltonian ${H}_{p}$ satisfies

$\Theta{H}_{p}(-t)\Theta=-{H}_{p}(t).$ (42)

Consider the formal expression of the evolution operator

$U(T,0)=\lim_{dt\rightarrow 0}\prod_{n=0}^{T/dt}e^{-i{H}_{p}(t_{n})dt},$

where $t_{n}=ndt$ and the product is time ordered. By inserting the resolution
of the identity $\hat{\mathbb{I}}=\Theta\Theta$ after each factor in the
product, and using relation (42), we get

$\displaystyle\Theta U(T,0)\Theta=$ $\displaystyle\lim_{dt\rightarrow
0}\prod_{n=0}^{T/dt}\Theta e^{-i{H}_{p}(t_{n})dt}\Theta$ $\displaystyle=$
$\displaystyle\lim_{dt\rightarrow
0}\prod_{n=0}^{T/dt}e^{i{H}_{p}(-t_{n})dt}=U^{\dagger}(0,-T)$

Let us take the trace of the evolution operator multiplied by $\Theta$, i.e.,

$\displaystyle\text{Tr}\left[U(T,-T)\Theta\right]=$
$\displaystyle\text{Tr}\left[U(0,-T)\Theta U(T,0)\right]$ $\displaystyle=$
$\displaystyle\text{Tr}\left[U(0,-T)\Theta U(T,0)\Theta\Theta\right]$
$\displaystyle=$
$\displaystyle\text{Tr}\left[U(0,-T)U^{\dagger}(0,-T)\Theta\right]=\text{Tr}\left[\Theta\right].$

This implies that the $S$ matrix satisfies
$\text{Tr}\left[S\Theta\right]=\text{Tr}\left[\Theta\right]=0$, which imposes
a constrain on the diagonal elements,

$S_{11}+S_{22}-S_{33}-S_{44}=0.$ (43)

(iv) Finally, the particular model (37) has an additional symmetry. By the
definition of the $S$ matrix, we have

$\displaystyle\hat{d}_{p}(+\infty)$
$\displaystyle=S_{31}\hat{c}_{p}+S_{32}\hat{c}_{-p}^{\dagger}+S_{33}\hat{d}_{p}+S_{34}\hat{d}_{-p}^{\dagger},$
$\displaystyle\hat{d}^{\dagger}_{-p}(+\infty)$
$\displaystyle=S_{41}\hat{c}_{p}+S_{42}\hat{c}_{-p}^{\dagger}+S_{43}\hat{d}_{p}+S_{44}\hat{d}_{-p}^{\dagger}.$

Taking the Hermitian conjugate of $\hat{d}_{p}(+\infty)$, and replacing $p$
with $-p$, we find

$\hat{d}^{\dagger}_{-p}(+\infty)=S^{*}_{31}(-p)\hat{c}^{\dagger}_{-p}+S^{*}_{32}(-p)\hat{c}_{p}+S^{*}_{33}(-p)\hat{d}_{-p}^{\dagger}+S^{*}_{34}(-p)\hat{d}_{p},$

where we used the fact that $S_{33}(-p)$ is real Li  _et al._ (2017).
Comparing the above two different expressions for $d^{\dagger}_{-p}(+\infty)$,
we find that

$S_{33}(-p)=S_{44}(p).$ (44)

Relations (40) - (44) are still insufficient for finding all elements of the
scattering matrix. Nevertheless, they allow us to find $\sum_{p}|S_{13}|^{2}$
without any approximation. Since the summation index in Eq. (18) in the main
text runs over a symmetric interval $p\in(-\pi,\pi)$, we have

$\displaystyle\sum_{p}|S_{13}|^{2}\stackrel{{\scriptstyle(\ref{one-
rel1})}}{{=}}$
$\displaystyle\sum_{p}\left[e^{-\frac{2\pi|\gamma|^{2}}{b_{1}+b_{2}}}-S_{11}S_{33}\right]$
(45) $\displaystyle\stackrel{{\scriptstyle(\ref{be2},\ref{one-rel4})}}{{=}}$
$\displaystyle\sum_{p>0}\left[2e^{-\frac{2\pi|\gamma|^{2}}{b_{1}+b_{2}}}-S_{11}(S_{33}+S_{44})\right]$
$\displaystyle\stackrel{{\scriptstyle(\ref{be2},\ref{one-rel3})}}{{=}}$
$\displaystyle
2\sum_{p>0}\left[e^{-\frac{2\pi|\gamma|^{2}}{b_{1}+b_{2}}}-|S_{11}|^{2}\right].$

##  Appendix C Bounds on $N^{o}_{\rm ex}$

In this section, we discuss the bounds on $N^{o}_{\rm ex}$, the number of
nonadiabatic excitations on the odd-site sublattice.

###  C.1 Lower bound

As given in the main text, the number of excitations on the odd-site
sublattice is expressed in terms of the $S$ matrix elements,

$N^{o}_{\rm ex}(+\infty)=\sum_{p}(|S_{33}|^{2}+|S_{31}|^{2}).$ (46)

The sum over the second term in the above equation can be obtained exactly
using the symmetry of the transition probability matrix (Eq. (13) in Ref. Li
_et al._ (2017)): $|S_{31}|=|S_{13}|$, which is given by Eq. (22) in the main
text.

The exact formula for $\sum_{p}|S_{33}|^{2}$ cannot be derived because energy
level $3$ does not have the highest slope. Nevertheless, a lower bound can be
found. Since $S_{33}$ and $S_{44}$ are both real,

$\displaystyle\sum_{p}|S_{33}|^{2}=$
$\displaystyle\sum_{p>0}(|S_{33}|^{2}+|S_{44}|^{2})$ (47) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{p>0}\left[(S_{33}+S_{44})^{2}+(S_{33}-S_{44})^{2}\right]$
$\displaystyle\geq$
$\displaystyle\frac{1}{2}\sum_{p>0}(S_{33}+S_{44})^{2}=\sum_{p}|S_{11}|^{2},$

where the last equality is a consequence of Eqs. (40) and (43). Thus,
$N^{o}_{\rm ex}$ is always bounded by the number of excitations on the even-
site sublattice, $N^{o}_{\rm ex}\geq N^{e}_{\rm ex}$.

###  C.2 Upper bound

The upper bound for the number of excitations of the odd-site sublattice can
be obtained under a special condition, namely, as a function of $p$, the sign
of $S_{33}$ remains fixed in a given regime of the quenching rate. In this
case, as a consequence of relation (44), $S_{33}S_{44}>0$. Hence,

$\sum_{p}|S_{33}|^{2}=\sum_{p>0}(|S_{33}|^{2}+|S_{44}|^{2})\leq\sum_{p>0}(S_{33}+S_{44})^{2}=\sum_{p>0}|2S_{11}|^{2}=2\sum_{p}|S_{11}|^{2},$
(48)

which is suppressed exponentially in the adiabatic limit. This condition is
fulfilled for large quenching rate, for which $S_{33}\sim 1$. We thus conclude
that for relatively large quenching rates, an exponential suppression for the
nonadiabatic excitation on the odd-site sublattice may be observed. However,
in the adiabatic limit, when the value of $S_{33}$ becomes small, this
condition may be broken, which makes it possible for the power-law excitation
tail to appear in the nearly adiabatic regime.

##  Appendix D Exact solutions for $J_{1}=J_{2}$

Finally, we note that the case with $J_{1}=J_{2}$ can be fully solved, i.e.,
all elements of the $S$ matrix can be obtained analytically. Let us group the
transition probabilities into a matrix with elements
$\hat{P}_{ij}\equiv|S_{ij}|^{2}$. For our 4$\times$4 model at $J_{1}=J_{2}$,
they are listed in Ref. Sinitsyn (2015):

$\hat{P}=\begin{pmatrix}p_{1}p_{2}&0&p_{2}q_{1}&q_{2}\\\
0&p_{1}p_{2}&q_{2}&p_{2}q_{1}\\\ p_{2}q_{1}&q_{2}&p_{1}p_{2}&0\\\
q_{2}&p_{2}q_{1}&0&p_{1}p_{2}\\\ \end{pmatrix},$ (49)

where

$p_{1}\equiv e^{-\frac{2\pi|g|^{2}}{|b_{1}-b_{2}|}},\ p_{2}\equiv
e^{-\frac{2\pi|\gamma|^{2}}{|b_{1}+b_{2}|}},\ q_{n}\equiv 1-p_{n},$ (50)

with restriction $b_{1}+b_{2}>0$, $b_{1}-b_{2}>0$. This solution can be
extended to a more general condition with $b_{1}\neq\pm b_{2}$, by permuting
the levels to full fill the condition. Now, with this exact solution, the
nonadiabatic excitations at $t=+\infty$ with the initial ground state
$|\downarrow\downarrow\downarrow\cdots\downarrow\downarrow\rangle$, Eqs. (18)
and (26) in the main text, can be computed as

$N^{e(o)}_{\rm
ex}(+\infty)=\sum_{p}e^{-\frac{2\pi|\gamma|^{2}}{|b_{1}+b_{2}|}},\
b_{1}\neq\pm b_{2},$ (51)

which confirms our analytical predictions in the main text for the case with
$J_{1}=J_{2}$. The above solution for the nonadiabatic excitations covers the
case of $b_{1}=b_{2}$, which corresponds to the simple uniform Ising chain.

##  Appendix E Nonadabatic excitations at $t=0$

The quasi-particle excitations at $t=0$ take the form of superpositions of
kinks. Here, we look at two different types of kinks separately, namely,
moving from left to right along the chain, the kinks between odd and even
sites, and the kinks between even and odd sites. The corresponding operators
are

$\displaystyle N^{e}_{\rm
kinks}=\frac{1}{2}\sum_{n}(1-\hat{\sigma}^{x}_{2n}\hat{\sigma}^{x}_{2n+1}),$
(52) $\displaystyle N^{o}_{\rm
kinks}=\frac{1}{2}\sum_{n}(1-\hat{\sigma}^{x}_{2n-1}\hat{\sigma}^{x}_{2n}).$

Using the Jordan-Wigner transformation (33) and the Fourier transformation
(34), the operators for the numbers of kinks can be expressed in terms of the
fermionic operators, i.e.,

$\displaystyle N^{e}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-(\hat{c}^{\dagger}_{p}\hat{d}_{p}+\hat{c}^{\dagger}_{p}\hat{d}^{\dagger}_{-p}+h.c.),$
(53) $\displaystyle N^{o}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-(\hat{c}^{\dagger}_{p}\hat{d}_{p}e^{-ip}+\hat{c}_{p}\hat{d}_{-p}e^{ip}+h.c.).$

Now the evolution of the fermionic operators in the Heisenberg picture is
governed by the effective Schrödinger equation with Hamiltonian $H_{p}$ in
(37). Define the $S$ matrix for the evolution over the time interval
$t\in(-T,0)$, i.e.,

$s=s(p)\equiv U^{p}(0,-T)_{T\rightarrow\infty},$ (54)

with $U^{p}(0,-T)$ the evolution matrix with the Hamiltonian $H_{p}$. The
fermionic operators in the Heisenberg picture at $t=0$ admit solutions in
terms of the $S$ matrix elements. We only list the relative quantities here:

$\displaystyle\hat{c}^{\dagger}_{p}(0)\hat{d}_{p}(0)$
$\displaystyle=s^{*}_{11}s_{31}\hat{c}^{\dagger}_{p}\hat{c}_{p}+s^{*}_{13}s_{33}\hat{d}^{\dagger}_{p}\hat{d}_{p}+\hat{O},$
(55) $\displaystyle\hat{c}^{\dagger}_{p}(0)\hat{d}^{\dagger}_{-p}(0)$
$\displaystyle=s^{*}_{11}s_{41}\hat{c}^{\dagger}_{p}\hat{c}_{p}+s^{*}_{13}s_{43}\hat{d}^{\dagger}_{p}\hat{d}_{p}+\hat{O},$
$\displaystyle\hat{c}_{p}(0)\hat{d}_{-p}(0)$
$\displaystyle=s_{12}s^{*}_{42}\hat{c}^{\dagger}_{-p}\hat{c}_{-p}+s_{14}s^{*}_{44}\hat{d}^{\dagger}_{-p}\hat{d}_{-p}+\hat{O}.$

Here, the operators on the RHS are interpreted as at $t=-\infty$ and the
operator $\hat{O}$ involves all other terms that do not contribute to the
computation of the number of kinks. Namely, we will evaluate the operator
expectation value with respect to the initial ground state at $t=-\infty$,
which is annihilated by $\hat{O}$. Taking average of the operators for the
number of kinks with respect to the state at $t=0$, we get the number of
kinks,

$\displaystyle N^{e}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-(s^{*}_{11}s_{31}+s^{*}_{13}s_{33}+s^{*}_{11}s_{41}+s^{*}_{13}s_{43}+c.g.),$
(56) $\displaystyle N^{o}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-(s^{*}_{11}s_{31}e^{-ip}+s^{*}_{13}s_{33}e^{-ip}+s_{12}s^{*}_{42}e^{ip}+s_{14}s^{*}_{44}e^{ip}+c.g.)$

Figure E1: The time-dependent spectrum of the Hamiltonian $H_{p}(t)$ from (37)
around the Dirac point (vertical red-dashed line), with the parameters
$p=0.3$, $b_{1}=2$, $b_{2}=1$, $J_{1}=3$, $J_{2}=1$.

The above involved $S$ matrix elements are represent in the diabatic basis,
i.e., the eigenbasis of Hamiltonian $H_{p}$ at $t=\infty$. One can express the
$S$ matrix element in terms of transition amplitude to the eigenbasis at
$t=0$, i.e., $\\{|\tilde{i}\rangle\\}_{i=1,2,3,4}$, as labeled in Fig. E1. For
instance,

$s_{11}=\langle 1|s|1\rangle=\langle
1|s\sum_{i}|\tilde{i}\rangle\langle\tilde{i}|1\rangle=\sum_{i}s_{1\tilde{i}}\langle\tilde{i}|1\rangle.$
(57)

Here the transition probability $s_{1\tilde{i}}$ can be solved using LZ
theory. Finally, in terms of transition amplitudes, the number of kinks reads

$\displaystyle N^{e}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-\frac{1}{4}\left[(e^{-ip}-1)(s_{3\tilde{3}}+s_{4\tilde{3}})-(e^{-ip}+1)(s_{3\tilde{4}}+s_{4\tilde{4}})+c.g.\right],$
(58) $\displaystyle N^{o}_{\rm
kinks}=\frac{1}{2}\sum_{p}1-\frac{1}{4}\left[(e^{-ip}-1)(e^{-ip}s_{3\tilde{3}}-e^{ip}s_{4\tilde{3}})-(e^{-ip}+1)(e^{-ip}s_{3\tilde{4}}-e^{ip}s_{4\tilde{4}})+c.g.\right].$

The excitation of kinks only involves transition amplitudes between the third
the fourth level around the Dirac point. This, as has been shown in the main
text for excitations at $t=+\infty$ on the odd-site sublattices, results in a
power-law scaling.


# A Measurement of CO(3–2) Line Emission from eBOSS Galaxies at $z\sim 0.5$
using Planck Data

Anirban Roy  0000-0001-5729-0246 Department of Physics, New York University,
726 Broadway, New York, NY, 10003, USA  Center for Computational Astrophysics,
Flatiron Institute, New York, NY 10010, USA  Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Nicholas Battaglia  0000-0001-5846-0411
Department of Astronomy, Cornell University, Ithaca, NY 14853, USA  Anthony R.
Pullen  Department of Physics, New York University, 726 Broadway, New York,
NY, 10003, USA  Center for Computational Astrophysics, Flatiron Institute, New
York, NY 10010, USA

###### Abstract

Line intensity mapping (LIM) is a novel observational technique in
astrophysics that utilizes the integrated emission from multiple atomic and
molecular transition lines from galaxies to probe the complex physics of
galaxy formation and evolution, as well as the large-scale structure of the
universe. Modeling multiple line luminosities of galaxies with varying masses
or their host halo masses poses significant uncertainty due to the lack of
observational data across a wide redshift range and the intricate nature of
astrophysical processes, making them challenging to model analytically or in
simulations. While future experiments aim to measure multiple line intensities
up to $z\sim 8$ across a wide volume using tomographic methods, we leverage
publicly available datasets from the CMB experiment Planck and the galaxy
survey eBOSS to constrain the CO(3-2) emission from galaxies. We correlate
galaxies from eBOSS data onto the full-sky CO(2-1) map produced by Planck and
report the first measurement of the average CO(3-2) intensity, $I_{CO}=45.7\pm
14.2\,\mathrm{Jy/sr}$ at $z\sim 0.5$ with $3.2\sigma$ confidence. Our findings
demonstrate that stacking methods are already viable with existing
observations from CMB experiments and galaxy surveys, and are complementary to
traditional LIM experiments.

Introduction.— Multi-probe cosmology integrates a diverse range of datasets to
study the properties and evolution of the universe, including those from
cosmic microwave background (CMB) radiation, large-scale structure (LSS)
surveys, supernovae, and galaxy clusters [1, 2]. Line intensity mapping (LIM)
emerges as a promising addition to this multifaceted approach, offering unique
insights into galaxy evolution and the large-scale distribution of matter in
the universe [3, 4, 5]. The ongoing and future LIM experiments aim to map
spatial intensity fluctuations of specific spectral lines, including
rotational level transitions of CO(1-0) to CO(13-12), [CII], [OIII],
H$\alpha$, H$\beta$, 21 cm and others, based on their brightness feature and
frequency overlaps [6, 7, 8, 9, 10, 11]. These observations across a broad
redshift range up to z $\sim 8$ could provide valuable understanding into the
formation of the first stars during the epoch of reionization, cosmic star
formation history, metal abundances in galaxies and cosmological information
such as the growth of structures, expansion history of the Universe and baryon
acoustic oscillations. While emissions of different J-level CO and [CII] from
individual galaxies are detected through high-resolution target observations
[12, 13, 14], the next frontier lies in detecting the spatial fluctuations of
these lines which is crucial for extracting both astrophysical and
cosmological information [15, 16].

One main challenge arises from modeling the multiple line intensities for
different types of galaxies across varying masses and environments, from a few
hundred million years after the Big Bang to the present day [17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Quantifying the various model
uncertainties is pivotal, as these could alter the interpretation of the
underlying physics of galaxy formation and evolution once the LIM observations
are conducted [32, 33]. Therefore, a promising approach is to consolidate all
probes related to galaxy formation and evolution that complement LIM into one
framework and perform a joint analysis to minimize these uncertainties.
Moreover, precise LIM observations are important for understanding foreground
contamination in CMB experiments, particularly for the measurement of
secondary anisotropies such as thermal and kinetic Sunyaev-Zeldovich effects
(tSZ and kSZ), lensing, and primordial $B$-mode signals, and this will aid in
mitigating such contamination for more accurate cosmological parameter
estimation [18, 34, 35].

Figure 1:  The panels depict intensity images generated by galaxies from four
distinct redshift bins, stacked onto the CO(2-1) Type 2 map derived from
Planck data. The panel corresponding to the redshift bin 0.49 to 0.57, aligned
with Planck’s frequency band of 219.3 to 230.8 GHz, is anticipated to exhibit
CO(3-2) emissions. Conversely, signals observed in the other panels reflect
the frequency-dependent residual foreground contamination. The FWHM of the
beam of the map is assumed to be $4^{\prime}.99$ [36].

The search for [CII] emission at $z\sim 2.6$ has been conducted through the
cross-correlation of Planck’s 353, 545, and 857 GHz data with BOSS quasars and
CMASS galaxies [37]. This analysis revealed an excess signal after jointly
fitting the [CII] intensity, CIB and tSZ model. The upper limit of CO(1-0)
emission from redshifts z = 2.4 to 3.4 has been estimated by cross-correlating
eBOSS quasars with the first-year intensity observations from COMAP [38]. In
this letter, we investigate the potential of employing a novel stacking method
to extract higher-order J-level CO transitions from the CO maps with lower
J-level. Drawing from the extensive use of stacking methods in detecting and
characterizing galaxy clusters on Planck’s tSZ maps [39, 40, 41, 42], we apply
this approach to existing component separated maps from Planck. Specifically,
we utilize the full-sky CO(2-1) map generated by Planck to identify the
contributions of CO(3-2) from galaxies detected by eBOSS at $z\sim 0.5$ [36,
43, 44]. The stacking of sources, such as galaxies in a map for our case
relies on minimizing the random noise properties while cumulatively adding the
weaker signals from each galaxy. This process increases the signal-to-noise
ratio (SNR) for the detection of the stacked image, enhancing its
detectability and then, inferring the average properties of sources in the
stacked image.

This letter is structured as follows: first, we introduce the methodology used
to correlate the galaxies observed by eBOSS with the Planck CO(2-1) map. Next,
we characterize potential contamination from other sources, such as residual
cosmic infrared background (CIB) and instrumental noise, and interpret the
stacked CO(3-2) intensity maps. Finally, we discuss the significance of these
results, highlighting the importance of these methods and their potential
implications for the current and next generation of LIM experiments. In this
letter, we adopt the cosmological model of a flat $\Lambda$CDM universe, as
established by the cosmological parameters derived from the Planck TT, TE,
EE+lowE+lensing results [45].

Cross-correlation Methodology.— Planck produced CO(1-0) and CO(2-1) maps at
the rest-frame frequencies of these lines, i.e., 115 GHz and 230 GHz, despite
lacking exact frequency channels to observe the sky at these frequencies [46].
This discrepancy was overcome due to the frequency mismatch of Planck’s
detectors, which extended their coverage to the rest-frame frequencies of
these two lines. Consequently, the Planck satellite was able to map these line
emissions with oversampling frequency bands of 109.6 – 115.4 GHz for CO(1-0)
and 219.3 – 230.8 GHz for CO(2-1) [36]. Line emissions from galaxies situated
at redshift $z$ to $z+\Delta z$, corresponding to the frequency bandwidth
$\nu$ to $\nu+\Delta\nu$, are expected to be redshifted to Planck’s detectors.
Additionally, higher-order CO transition lines from galaxies at earlier cosmic
times are also redshifted to the same frequency channels. Therefore, Planck CO
maps not only comprise prominent CO(1-0) and CO(2-1) lines but also include
the subdominant contributions from higher-order CO lines. For instance,
CO(3-2) emissions from redshift 0.5724 to 0.4940 are redshifted to Planck’s
oversampling frequency 219.3 – 230 GHz, and CO(2-1) from redshift 1.1 to 0.99
are redshifted to Planck’s CO(1-0) map. With this premise, we decide to look
for the CO(3-2) signal from Planck’s CO(2-1) maps.

Planck data yield two distinct map types: Type 1 and Type 2 [36]. Type 1 maps,
derived from spectral transitions of bolometer pairs within a single frequency
channel, do not capture intensity contamination from other frequency channels.
In contrast, Type 2 maps integrate intensity responses across multiple
frequency channels, facilitating the separation of foreground emissions like
CMB, free-free, and dust emissions from CO emissions. Comparison of the
standard deviations at high latitudes reveals that the uncertainty in Type 2
maps is approximately four times smaller than that in Type 1 maps [36].
Consequently, unless specified otherwise, we use CO(2-1) Type 2 maps for our
stacking analysis aimed at detecting CO(3-2) signals at redshift $z\sim 0.5$.

Within the redshift range from $z=0.57-0.49$, we utilize the galaxy catalogs
of eBOSS to ensure alignment of CO(3-2) emission from these galaxies with the
exact frequency bandwidth used in constructing Planck’s CO(2-1) map. A total
of $427,323$ galaxies within this redshift bin are identified for stacking
analysis. Using the Healpy package [47], we identify the positions of
individual galaxies on the CO(2-1) map and extract circular regions with a
$1.5\degree$ radius, centered on each galaxy’s right ascension (RA) and
declination (DEC) coordinates. We select this radius to ensure that the
resulting map remains within a square area enclosed by the circle, with a side
length of $2.12\degree$. We divided this square region into 125 pixels,
ensuring that each pixel corresponds to an approximate size of $\sim
1^{\prime}$. We iterate this process for all the galaxies at this redshift bin
and finally, we stack all those cutout maps and make a single stacked image
that captures the signature of CO(3-2) emissions from all the galaxies.

During the stacking analysis, the presence of any correlated signal, rather
than random noise, is anticipated to appear at the centre of the image, given
that emissions primarily originate from galaxies. However, such observations
do not definitively ascertain whether the signal originates solely from
CO(3-2) emission or from residual foregrounds correlated with the galaxy
positions. To explore this, we employ a stacking approach on galaxies from
adjacent redshift bins. This allows us to ascertain if CO(3-2) emission is
absent in stacked images composed of galaxies from other redshifts. To ensure
consistency in our analysis, an equal number of galaxies are selected from
subsequent redshift bins from which we do not expect CO(3-2) emission.

Subsequently, we identify adjacent redshift bins to $z_{i}$ and $z_{i}+\Delta
z_{i}$ in a manner ensuring that each bin contains the same number of
galaxies. Similar to the redshift bin of CO(3-2) as described previously, we
conduct stacking analysis on these adjacent bins using the Planck CO(2-1) Type
2 map.

We present the result of the stacking analysis at four redshifts in Figure
@stylefig:stacked_maps. We do not find any significant signal in the stacked
image of galaxies situated between redshifts 0.57 and 0.70, as expected due to
the absence of CO(3-2) emission at these redshifts. The emissions in the
redshift range of 0.49 to 0.57 indicate where CO(3-2) is expected to appear,
shown in the second panel from the left. Furthermore, as we reach to the low
redshift bins, we do not see the signal disappear which would be expected if
the intensity signal came purely from CO(3-2) emission. In contrast, we see
the mean amplitude of the signal at the center of the stacked map increases
significantly, confirming galaxies present in the maps, and CO(3-2) signal
adds up with the residual signal only on the redshift bin 0.49 to 0.57. It is
evident that the size of the central region increases due to the fact that the
angular size of galaxies increases at the low redshifts than the high
redshift, hence the intensity distribution spreads over more pixels at low
redshift than at high redshifts.

We assess the intensity of galaxies by employing an aperture photometry
technique facilitated by the Astropy package [48]. Using a central circle
radius of 10 arcminutes, we capture the primary intensity concentration
attributed to the stacked galaxies. Initially, we determine the intensity of
this concentration and subsequently compute the variance of the intensity
measurement relative to the background intensity distribution. To accomplish
this, we randomly select 100 circles outside the central zone (10 arcminutes)
and compute the intensity variance of these outer regions. This approach
allows us to quantify both the intensity of emissions originating from
galaxies and the errors stemming from random background noise.

We conduct several tests to ensure the efficacy of the stacking analysis on
our target galaxies and to ascertain the appropriateness of signals from
various redshift bins for a model encompassing CO(3-2) signal and residual
foregrounds. For the initial test, we use the same number of galaxies
corresponding to the redshift bin for the emission of CO(3-2) and randomly
disperse them across the sky. Due to variations in their RA and DEC values,
galaxies in the eBOSS catalog do not match this distribution. Afterwards, we
perform the stacking analysis on these randomly distributed galaxies projected
to Planck’s CO(2-1) Type 2 map as described previously. Our results reveal a
lack of observable structures at the center of the stacked map, but trace the
presence of homogeneous random noise across the entire map. Moreover, the
overall amplitude is approximately three times smaller than the intensity
observed at the core of the stacked eBOSS galaxies on the CO(2-1) map. This
outcome confirms that the signal observed in Figure @stylefig:stacked_maps
indeed originates from the eBOSS galaxies within the selected redshift bins.

Figure 2:  We display the estimated intensities along with their respective
errors with boxes at four redshift bins. The red line represents the residual
CIB spectrum calculated using the best-fit parameters. The shaded green region
indicates the redshift range where we anticipate contributions from CO(3-2),
while the shaded red regions indicate the residual CIB emissions at those
redshift bins.

Interpretation.— We proceed to quantify the remaining foreground in the
stacked maps of CO(3-2) across four distinct redshift bins, in order to
understand the emission we are seeing in these bins. Calculating the
intensities of the stacked image across four adjacent redshift bins provides
insights into how residual CIB varies with observed frequency (or redshift).
It is anticipated that the CIB spectrum should closely resemble a blackbody
spectrum; however, the amplitude is expected to be reduced due to the
component separation process [36]. We model the residual CIB emission using a
power-law parameterization as $I_{\rm
res}(z)=I^{0}_{res}\times(1+z)^{-\alpha}$, where $I^{0}_{res}$ is the
amplitude of residual foreground at $z\sim 0$ and $\alpha$ is the power-law
index.

Figure 3: The constraints on parameters for residual CIB and average CO(3-2)
luminosity at $z\sim 0.5$ estimated from the stacking analysis with eBOSS
galaxies onto the Planck’s CO(2-1) map.

To ascertain the validity of this assumption for the foreground modelling, we
performed stacking analysis using an equal number of galaxies across four
redshift bins with a CO(1-0) map, as opposed to using CO(2-1) map as
previously described. Due to a frequency mismatch, we do not expect any CO
emissions for the same redshift bins to be present in the CO(1-0) map. The
results, detailed in Figure @stylefig:sub_CIB in the supplementary section,
reveal a diminished signal, with no prominent structures discernible at the
center of the stacked image. However, upon fitting the intensity data points
for this scenario using a model considering only residual CIB, we derived
$I^{0}_{res}=133^{+14}_{-12}\,(Jy/sr)$ and $\alpha=0.52^{+0.25}_{-0.24}$. In
this instance, the chi-square value is 0.55 with 2 degrees of freedom,
yielding a probability-to-exceed (PTE) value of 0.76. As this model can
perfectly fit the residual foregrounds at the 115 GHz map, we apply the same
foreground model to probe the CO(3-2) signal in the 217 GHz channel.

Parameter |  without CO |  with CO
---|---|---
$I_{0}$ |  $545^{+56}_{-51}$ |  $587^{+63}_{-59}$
$\alpha$ |  $1.94^{+0.29}_{-0.27}$ |  $2.22\pm 0.29$
$L_{\rm CO(3-2)}$ |  N/A |  $1495^{+460}_{-465}$
$\chi^{2}$ |  7.95 |  1.31
PTE |  0.018 |  0.25
Table 1: Best fit parameters with 68% confidence intervals, accompanied by the
corresponding values of $\chi^{2}$ and PTE, obtained from fitting the
intensities of stacked profiles at four redshifts. The fits consider
contributions from residual CIB alone and residual CIB with CO(3-2) emissions.

We construct a model that combines both residual CIB and CO(3-2) emissions,
expressed as $I_{\rm stack}=I_{\rm res}+I_{\rm CO(3-2)}$. The CO(3-2) emission
model is assumed to originate from the specific redshift bin of 0.49-0.57;
therefore, we introduce a free parameter $L_{\rm CO(3-2)}$, representing the
average CO(3-2) luminosity of galaxies within that bin. This parameter relates
to the intensity as $I_{\rm CO(3-2)}=L_{\rm CO(3-2)}\times(c/4\pi\nu_{\rm
rest}H(z))$, where $\nu_{\rm rest}$ is the rest-frame frequency of CO(3-2) and
$H(z)$ is the Hubble parameter at the redshift of CO(3-2) emission. To measure
the average intensity of CO(3-2), we integrate $dn/dz$ of galaxies in redshift
bins and compute the average intensity or foreground as $\langle
I\rangle=\int_{z_{\rm min}}^{z_{\rm max}}I(z)\frac{dn(z)}{dz}dz/\int_{z_{\rm
min}}^{z_{\rm max}}\frac{dn}{dz}dz$. Here, $I(z)$ is the total intensity at a
given redshift, including residual foreground and CO(3-2) signal, and $dn/dz$
is the redshift distribution of eBOSS galaxies.

In Figure @stylefig:co32-detection, we illustrate how the intensity of the
stacked images, encompassing CO(3-2) emissions around $z\sim 0.5$ and residual
CIB, varies with redshift. As expected, the intensity of the residual CIB
decreases at the high redshift showing a power-law trend. We estimate the
intensity at the redshift bins 0.18–0.37, 0.37–0.49, and 0.57–0.70 to be
$339.6\pm 15.1$, $270\pm 9.2$, and $195.5\pm 12$ Jy/sr, and the intensity in
the redshift bin where CO(3-2) contributes to be $266\pm 11.5$ Jy/sr. We make
use of the emcee package [49] to fit the parameters of the intensity model to
the estimated intensities at four redshift bins. We consider two scenarios: i)
fitting the data while considering only residual CIB present in the stacked
map, where the parameters are $I^{0}_{res}$ and $\alpha$, and ii) considering
both residual CIB and CO at redshift bin $z\sim 0.5$, where we introduce an
additional parameter $L_{\rm CO(3-2)}$ and subsequently calculate the
intensity of CO(3-2) from luminosity. We show the constraints on parameters
through a triangle plot in Figure @stylefig:triangle_plot and the best-fit
parameters for both cases are reported in Table
@styletab:parameter_comparison. When fitting the data with a CIB-only model,
the $\chi^{2}$ becomes 7.95 with 2 degrees of freedom, yielding a PTE value of
0.018. Conversely, fitting the data with residual CIB and CO(3-2) emission
reduces the $\chi^{2}$ value to 1.31, with a corresponding PTE value of 0.25.
Therefore, our analysis confirms that the presence of CO(3-2) emission, along
with residual foreground, provides a better description of the observed data.

Figure 4: Various CO(3-2) model predictions are compared with the measured
CO(3-2) intensity at $z\sim 0.5$. The legend provides descriptions of the SFR
models, while annotations within the gray boxes denote the models that convert
SFR to CO(3-2) line luminosities. The horizontal bar for a particular line
luminosity model is calculated from the variance in intensities for the
different SFR models.

This measurement translates to an average CO(3-2) luminosity density of
$L_{\rm CO}=1495^{+460}_{-465}\,L_{\odot}/{Mpc}^{3}$ at $z\sim 0.5$, with the
corresponding intensity constrained to $I_{\rm CO}=45.7\pm 14.2$ with 68%
confidence. Predicting the intensity of CO(3-2) analytically is challenging
due to various factors: i) the relationship between CO(3-2) luminosity and
star formation rate (SFR) or halo mass is not precisely constrained and
depends on multiple empirical scaling relations, ii) the ratio of CO(3-2) to
CO(1-0) or other CO J-level transitions is not firmly established across a
wide redshift range, iii) observed CO(3-2) intensity may deviate from
estimated intensity due to differences in integrating luminosity over all halo
masses versus specific galaxy selection criteria, among other factors.
Considering these complexities, we compare our measurement with the estimated
values with those in the literature using different SFR models and scaling
relations between SFR and CO(3-2) luminosity, Visbal10 [50], Greve14 [12], and
Kamenetzky15 [51]. We estimate the average intensity of CO(3-2) lines for
eBOSS galaxies using the formula
$I_{\text{est}}=\frac{c}{4\pi\nu_{\text{rest}}H(z=0.5)}\times\int_{M_{\text{min}}}^{M_{\text{max}}}L(M)\frac{dn}{dM}dM$.
Here, $M_{\text{min}}$ and $M_{\text{max}}$ represent the minimum and maximum
stellar masses of the galaxy samples, respectively. We present our findings in
Figure @stylefig:model_comparison. Details of the model description and
theoretical background are provided in [33]. We calculate the average
intensity of CO(3-2) using five SFR models: Tng100, Tng300 [52, 53, 54],
Fonseca16 [20], Silva15 [55], and Behroozi19 [56] as described in [33]. The
variance and mean of these intensities corresponding to SFR models are
determined, and the error is calculated from the variance. Our measurement is
in good agreement with Visbal10 and Greve14, as most of these models fall
within the $1\sigma$ uncertainty of our measurement.

Conclusion and Outlook.— In this letter, we investigate the feasibility of
employing stacking analysis as a novel method to constrain the average
intensity of line emissions. We present the first constraints on the
luminosity and intensity of CO(3-2) by employing this method on existing data
from Planck and eBOSS. We constrained $I_{CO}=45.7\pm 14.2\,\mathrm{Jy/sr}$ at
$z\sim 0.5$ with $3.2\sigma$ confidence. Comparing our results to three
empirical models of CO(3-2) emission we found that they were in agreement with
our measurement, but this agreement depended upon the halo mass to star-
formation relation chosen. When used on forthcoming line intensity mapping
datasets, we expect a higher SNR yield for similar cross-correlation analyses,
and such analyses will have the potential to further our understanding of
galaxy evolution [57, 38, 58]. Forthcoming LIM cross-correlation analyses like
this will be able to narrow down the predictions for extragalactic CO
contamination in broadband CMB observations [e.g., 35, 59].

The method discussed in this letter complements the more commonly discussed
LIM analysis at the power spectrum level. Analyzing LIM at the power spectrum
has the potential to constrain the mean intensity of lines in a tomographic
manner, aiding in constraining the average line intensity. However, a
degeneracy exists between the average intensity and bias parameters, with
ongoing experiments able to place limits on the combined of two, $b\times
I_{\nu}$ parameter. In contrast, $I_{\nu}$ can be estimated solely through
stacking analysis for a given sample of galaxies, independent of the bias
parameter. Therefore, assimilating both probes—stacking and power spectrum
analysis—could lead to tighter constraints on the average intensity parameter
as well as other astrophysical parameters.

Our work is a first step towards detecting the intensity of a particular line
by correlating galaxies with an intensity map and developing the pipeline for
data analysis. In the future, an extensive array of galaxy samples will become
available from experiments such as DESI, LSST, and Euclid, accompanied by low-
noise LIM maps from experiments like EoR-Spec at FYST, COMAP, and CONCERTO,
thereby expanding the scientific potential of stacking analysis. Following the
reported detection of the stacking of CO(3-2) lines of galaxies in this
letter, we plan to present a detailed analysis of robust forecasts for future
experiments in our subsequent works.

We thank Dongwoo Chung, Colin Hill, Abhishek Maniyar, Matt McQuinn, and Alex
van Engelen for their helpful discussions. AR acknowledges support from NASA
under award number 80NSSC18K1014939 during the final stage of this work. NB
acknowledges support from CCAT collaboration for this work and additional
support from NASA grants 80NSSC18K0695 and 80NSSC22K0410. ARP was supported by
NASA under award numbers 80NSSC18K1014, NNH17ZDA001N, and 80NSSC22K0666, and
by the NSF under award number 2108411.

## References

  * Abbott et al. [2019] T. Abbott, A. Alarcon, S. Allam, P. Andersen, F. Andrade-Oliveira, J. Annis, J. Asorey, S. Avila, D. Bacon, N. Banik, et al., Physical review letters 122, 171301 (2019).
  * Moresco et al. [2022] M. Moresco, L. Amati, L. Amendola, S. Birrer, J. P. Blakeslee, M. Cantiello, A. Cimatti, J. Darling, M. Della Valle, M. Fishbach, et al., Living Reviews in Relativity 25, 6 (2022), eprint 2201.07241.
  * Kovetz et al. [2017] E. D. Kovetz, M. P. Viero, A. Lidz, L. Newburgh, M. Rahman, E. Switzer, M. Kamionkowski, J. Aguirre, M. Alvarez, J. Bock, et al., arXiv preprint arXiv:1709.09066 (2017).
  * Kovetz et al. [2020] E. D. Kovetz et al., Bull. Am. Astron. Soc. 51, 101 (2020), eprint 1903.04496.
  * Schaan and White [2021] E. Schaan and M. White, JCAP 05, 068 (2021), eprint 2103.01964.
  * Crites et al. [2014] A. T. Crites, J. J. Bock, C. M. Bradford, T. C. Chang, A. R. Cooray, L. Duband, Y. Gong, S. Hailey-Dunsheath, J. Hunacek, P. M. Koch, et al., in _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII_ , edited by W. S. Holland and J. Zmuidzinas (2014), vol. 9153 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 91531W.
  * Doré et al. [2018] O. Doré, M. W. Werner, M. L. Ashby, L. E. Bleem, J. Bock, J. Burt, P. Capak, T.-C. Chang, J. Chaves-Montero, C. H. Chen, et al., arXiv preprint arXiv:1805.05489 (2018).
  * Ade et al. [2020a] P. Ade et al. (CONCERTO), Astron. Astrophys. 642, A60 (2020a), eprint 2007.14246.
  * Ade et al. [2020b] P. A. R. Ade et al., J. Low Temp. Phys. 199, 1027 (2020b), eprint 1912.07118.
  * Aravena et al. [2022] M. Aravena, J. E. Austermann, K. Basu, N. Battaglia, B. Beringue, F. Bertoldi, F. Bigiel, J. R. Bond, P. C. Breysse, C. Broughton, et al., The Astrophysical Journal Supplement Series 264, 7 (2022).
  * Cleary et al. [2022] K. A. Cleary, J. Borowska, P. C. Breysse, M. Catha, D. T. Chung, S. E. Church, C. Dickinson, H. K. Eriksen, M. K. Foss, J. O. Gundersen, et al., The Astrophysical Journal 933, 182 (2022).
  * Greve et al. [2014] T. Greve, I. Leonidaki, E. Xilouris, A. Weiss, Z.-Y. Zhang, P. van der Werf, S. Aalto, L. Armus, T. Diaz-Santos, A. Evans, et al., arXiv preprint arXiv:1407.4400 (2014).
  * Genzel et al. [2015] R. Genzel, L. Tacconi, D. Lutz, A. Saintonge, S. Berta, B. Magnelli, F. Combes, S. García-Burillo, R. Neri, A. Bolatto, et al., The Astrophysical Journal 800, 20 (2015).
  * Schaerer et al. [2020] D. Schaerer, M. Ginolfi, M. Béthermin, Y. Fudamoto, P. Oesch, O. Le Fèvre, A. Faisst, P. Capak, P. Cassata, J. Silverman, et al., Astronomy & Astrophysics 643, A3 (2020).
  * Silva et al. [2021] M. B. Silva, E. D. Kovetz, G. K. Keating, A. M. Dizgah, M. Bethermin, P. C. Breysse, K. Karkare, J. L. Bernal, and J. Delabrouille, Experimental Astronomy 51, 1593 (2021).
  * Padmanabhan [2021] H. Padmanabhan, International Journal of Modern Physics D 30, 2130009 (2021).
  * Suginohara et al. [1998] M. Suginohara, T. Suginohara, and D. Spergel, The Astrophysical Journal 512, 547 (1998).
  * Righi et al. [2008] M. Righi, C. Hernandez-Monteagudo, and R. Sunyaev, Astron. Astrophys. 489, 489 (2008), eprint 0805.2174.
  * Carilli [2011] C. Carilli, The Astrophysical Journal Letters 730, L30 (2011).
  * Fonseca et al. [2017] J. Fonseca, M. B. Silva, M. G. Santos, and A. Cooray, Mon. Not. Roy. Astron. Soc. 464, 1948 (2017).
  * Gong et al. [2017] Y. Gong, A. Cooray, M. B. Silva, M. Zemcov, C. Feng, M. G. Santos, O. Dore, and X. Chen, The Astrophysical Journal 835, 273 (2017).
  * Chung et al. [2020] D. T. Chung, M. P. Viero, S. E. Church, and R. H. Wechsler, Astrophys. J. 892, 51 (2020), eprint 1812.08135.
  * Padmanabhan [2018] H. Padmanabhan, Mon. Not. Roy. Astron. Soc. 475, 1477 (2018), eprint 1706.01471.
  * Padmanabhan [2019] H. Padmanabhan, Mon. Not. Roy. Astron. Soc. 488, 3014 (2019), eprint 1811.01968.
  * Dumitru et al. [2019] S. Dumitru, G. Kulkarni, G. Lagache, and M. G. Haehnelt, Mon. Not. Roy. Astron. Soc. 485, 3486 (2019), eprint 1802.04804.
  * Chung et al. [2019] D. T. Chung et al. (COMAP), Astrophys. J. 872, 186 (2019), eprint 1809.04550.
  * Kannan et al. [2022] R. Kannan, A. Smith, E. Garaldi, X. Shen, M. Vogelsberger, R. Pakmor, V. Springel, and L. Hernquist, Mon. Not. Roy. Astron. Soc. 514, 3857 (2022), eprint 2111.02411.
  * Murmu et al. [2021] C. S. Murmu, S. Majumdar, and K. K. Datta, Monthly Notices of the Royal Astronomical Society 507, 2500 (2021).
  * Karoumpis et al. [2022] C. Karoumpis, B. Magnelli, E. Romano-Díaz, M. Haslbauer, and F. Bertoldi, Astronomy & Astrophysics 659, A12 (2022).
  * Sun et al. [2023] G. Sun, L. Mas-Ribas, T.-C. Chang, S. R. Furlanetto, R. H. Mebane, M. O. Gonzalez, J. Parsons, and A. Trapp, The Astrophysical Journal 950, 40 (2023).
  * Garcia et al. [2023] K. Garcia, D. Narayanan, G. Popping, R. Anirudh, S. Sutherland, and M. Kaasinen, arXiv preprint arXiv:2311.01508 (2023).
  * Murmu et al. [2023] C. S. Murmu, K. P. Olsen, T. R. Greve, S. Majumdar, K. K. Datta, B. R. Scott, T. D. Leung, R. Davé, G. Popping, R. O. Ochoa, et al., Monthly Notices of the Royal Astronomical Society 518, 3074 (2023).
  * Roy et al. [2023] A. Roy, D. Valentín-Martínez, K. Wang, N. Battaglia, and A. van Engelen, arXiv preprint arXiv:2304.06748 (2023).
  * Puglisi et al. [2017] G. Puglisi, G. Fabbian, and C. Baccigalupi, Mon. Not. Roy. Astron. Soc. 469, 2982 (2017), eprint 1701.07856.
  * Maniyar et al. [2023] A. S. Maniyar, A. Gkogkou, W. R. Coulton, Z. Li, G. Lagache, and A. R. Pullen, Physical Review D 107, 123504 (2023).
  * Ade et al. [2014a] P. A. Ade, N. Aghanim, M. Alves, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. Banday, et al., Astronomy & Astrophysics 571, A13 (2014a).
  * Pullen et al. [2018] A. R. Pullen, P. Serra, T.-C. Chang, O. Dore, and S. Ho, Mon. Not. Roy. Astron. Soc. 478, 1911 (2018), eprint 1707.06172.
  * Dunne et al. [2024] D. A. Dunne, K. A. Cleary, P. C. Breysse, D. T. Chung, H. T. Ihle, J. R. Bond, H. K. Eriksen, J. O. Gundersen, L. C. Keating, J. Kim, et al., The Astrophysical Journal 965, 7 (2024).
  * Greco et al. [2015] J. P. Greco, J. C. Hill, D. N. Spergel, and N. Battaglia, The Astrophysical Journal 808, 151 (2015).
  * Aghanim et al. [2016] N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. Banday, R. Barreiro, J. Bartlett, N. Bartolo, E. Battaner, et al., Astronomy & Astrophysics 594, A22 (2016).
  * Hill et al. [2018] J. C. Hill, E. J. Baxter, A. Lidz, J. P. Greco, and B. Jain, Physical Review D 97, 083501 (2018).
  * Hall et al. [2019] K. R. Hall, N. L. Zakamska, G. E. Addison, N. Battaglia, D. Crichton, M. Devlin, J. Dunkley, M. Gralla, J. C. Hill, M. Hilton, et al., Monthly Notices of the Royal Astronomical Society 490, 2315 (2019).
  * Dawson et al. [2016] K. S. Dawson, J.-P. Kneib, W. J. Percival, S. Alam, F. D. Albareti, S. F. Anderson, E. Armengaud, É. Aubourg, S. Bailey, J. E. Bautista, et al., The Astronomical Journal 151, 44 (2016).
  * Blanton et al. [2017] M. R. Blanton, M. A. Bershady, B. Abolfathi, F. D. Albareti, C. A. Prieto, A. Almeida, J. Alonso-García, F. Anders, S. F. Anderson, B. Andrews, et al., The Astronomical Journal 154, 28 (2017).
  * Aghanim et al. [2020] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. Barreiro, N. Bartolo, S. Basak, et al., Astronomy & Astrophysics 641, A6 (2020).
  * Ade et al. [2014b] P. A. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. Banday, R. Barreiro, et al., Astronomy & Astrophysics 571, A9 (2014b).
  * Gorski et al. [2005] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, The Astrophysical Journal 622, 759 (2005).
  * Price-Whelan et al. [2022] A. M. Price-Whelan, P. L. Lim, N. Earl, N. Starkman, L. Bradley, D. L. Shupe, A. A. Patil, L. Corrales, C. Brasseur, M. Nöthe, et al., The Astrophysical Journal 935, 167 (2022).
  * Foreman-Mackey et al. [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Visbal and Loeb [2010] E. Visbal and A. Loeb, Journal of Cosmology and Astroparticle Physics 2010, 016 (2010).
  * Kamenetzky et al. [2016] J. Kamenetzky, N. Rangwala, J. Glenn, P. Maloney, and A. Conley, The Astrophysical Journal 829, 93 (2016).
  * Nelson et al. [2018] D. Nelson et al. (2018), eprint 1812.05609.
  * Pillepich et al. [2018] A. Pillepich et al., Mon. Not. Roy. Astron. Soc. 475, 648 (2018), eprint 1707.03406.
  * Springel et al. [2018] V. Springel et al., Mon. Not. Roy. Astron. Soc. 475, 676 (2018), eprint 1707.03397.
  * Silva et al. [2015] M. B. Silva, M. G. Santos, A. Cooray, and Y. Gong, Astrophys. J. 806, 209 (2015), eprint 1410.4808.
  * Behroozi et al. [2019] P. Behroozi, R. H. Wechsler, A. P. Hearin, and C. Conroy, Mon. Not. Roy. Astron. Soc. 488, 3143 (2019), eprint 1806.07893.
  * Breysse et al. [2022] P. C. Breysse, S. Yang, R. S. Somerville, A. R. Pullen, G. Popping, and A. S. Maniyar, The Astrophysical Journal 929, 30 (2022).
  * Horlaville et al. [2024] P. Horlaville, D. T. Chung, J. R. Bond, and L. Liang, Monthly Notices of the Royal Astronomical Society 531, 2958 (2024).
  * Kokron et al. [2024] N. Kokron, J. L. Bernal, and J. Dunkley, arXiv e-prints arXiv:2405.20369 (2024), eprint 2405.20369.

##  Appendix A Characterization of residual foreground

We present a methodology for confirming the correlation between CO(3-2)
emissions and galaxies within the frequency bandwidth observed by Planck.
Specifically, we performed checks and tests to ensure that the CO(3-2) signal
originates from galaxies corresponding to Planck’s frequency range in the
presence of residual foreground on the map. We emphasize that within the
observational frequency range of 109.6–115.4 GHz, where the Planck CO(1-0) map
has been constructed, there should not be correlated CO(3-2) emissions from
galaxies in the redshift bin of 0.49 to 0.57. To verify this, we employ a
stacking technique across galaxies grouped into four distinct redshift bins,
aligning with the Planck CO(1-0) map. This approach enables us to effectively
characterize and account for residual foreground contamination present in the
map at 115 GHz.

We show the stacking images at four redshifts in Figure @stylefig:sub_CIB,
encompassing residual foreground emissions at 115 GHz. As CO emission for all
J-level transitions should not be present at any of these four redshifts, we
expect these images to capture the information of the residual foreground. As
previously mentioned, we modeled the residual foreground as a power-law with
two free parameters, $I^{0}_{\rm res}$ and $\alpha$. We aim to check the
validity of this model by estimating the $\chi^{2}$ and PTE for the 115 GHz
channel. The average intensity changes from $116.2\pm 5.3$ Jy/sr to $101.3\pm
4.6$ Jy/sr from the redshift bins 0.18–0.37 to 0.57–0.70. We show the
distribution of mean intensity with redshift in Figure @stylefig:CIB-res-fit-
CO10, which displays a power-law trend. We estimate the best-fit values of
$I_{\rm res}^{0}$ and $\alpha$ to be $133^{+14}_{-12}\,(Jy/sr)$ and
$0.52^{+0.25}_{-0.24}$, respectively. The $\chi^{2}$ value is 0.55 with 2
degrees of freedom, and the PTE is 0.76. This test allows us to assume a
power-law-like residual foreground model for the CO(2-1) map constructed at
230 GHz to probe the CO(3-2) signal from the galaxies at $z\sim 0.5$.

Figure 5:  The left-over contaminants in the map within the observational
frequency bandwidth of 109.6 – 115.4 GHz. The absence of prominent structure
at the center of the stacked images signifies that there is no CO(3-2)
contribution from galaxies in the redshift bin 0.49 to 0.57, and also that the
residual CIB contributions are smaller than those for the observational
frequency range of 219.3 – 230.8 GHz. The assumed FWHM of the map’s beam is
$9^{\prime}.65$ [36]. Figure 6: Mean intensities of stacked eBOSS galaxies
overlaid onto Planck’s CO(1-0) map. The boxes represent intensity variations,
with box length indicating intensity uncertainty, and box width depicting the
range of redshift bins. Additionally, the figure displays the best-fit power-
law model for the residual foreground, aligning with the measured intensities.