††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Baryogenesis via QCD preheating with nonadiabatic baryon chemical potential

Jimin Wang  Center for Theoretical Physics and College of Physics, Jilin
University, Changchun, 130012, China.  Xin-Ru Wang  Center for Theoretical
Physics and College of Physics, Jilin University, Changchun, 130012, China.
Shinya Matsuzaki  Center for Theoretical Physics and College of Physics, Jilin
University, Changchun, 130012, China.

###### Abstract

The chiral phase transition in QCD can be supercooled in the thermal history
of the universe to be instantaneously out-of equilibrium, if QCD is coupled to
a dark QCD sector exhibiting the dark chiral phase transition of the first
order. In that case the QCD sigma meson field (as the chiral order parameter,
or the light quark condensate) starts to roll in a nonadiabatic way down to
the true QCD vacuum. Meanwhile a dynamic baryonic chemical potential can be
generated solely within QCD, which is governed by the dynamic motion of the
QCD sigma meson field, analogously to the spontaneous baryogenesis or the
leptogenesis via the Higgs or axionlike relaxation scenario. When QCD is
further allowed to communicate with a dark fermion with mass of order of 1 GeV
and the baryon number violating coupling to neutron, the nonadiabatic QCD
sigma motion along with the nonadiabatic baryon chemical potential can trigger
the preheating and produce the baryon number asymmetry. We discuss this
scenario in details to find that the QCD-induced dynamic baryon chemical
potential plays a significant role for the QCD preheating and the
baryogenesis, which yields the desired amount of the asymmetry today
consistently with current astrophysical, cosmological, and terrestrial
experimental constraints. Cosmological and phenomenological consequences
characteristic to the present scenario are also addressed.

##  I Introduction

The cosmic microwave background (CMB) and the big bang nucleosynthesis (BBN)
imply existence of the baryon number asymmetry (BNA) today, which is captured
in terms of the ratio of the baryon number density to the entropy density with
the amount $(n_{B}-\bar{n}_{B})/s\sim 6\times 10^{-10}$ Planck:2015fie ;
Cyburt:2015mya ; Planck:2018vyg , indicating the dominance of matters against
anti-matters. The Standard Model (SM) of particle physics cannot sufficiently
produce the present-day BNA, so this BNA issue is one of well motivated
slogans for theoretical particle physicists to call for beyond the SM.

The mechanism to generate this asymmetry, called the baryogenesis, needs to
satisfy the Sakharov’s criteria Sakharov:1967dj : C and CP violation, baryon
number violation, and out-of-equilibrium. Among the Sakharov’s criteria, in
particular, realization of non-equilibrium can involve a cosmological aspect,
i.e., a dynamic and nonadiabatic time evolution of the vacuum. Such a dynamic
aspect of the vacuum has been discussed aiming at the post inflationary
universe, so-called preheating Dolgov:1989us ; Traschen:1990sw ; Kofman:1994rk
; Shtanov:1994ce ; Kofman:1997yn  (for reviews, see e.g., Kofman:1997yn ;
Amin:2014eta ; Lozanov:2019jxc ). When the inflaton couples to the SM
particles, the nonadiabatic motion or oscillation of the inflaton generates
the nonadiabatic quantum states for the SM particles, that is what is called
the nonperturbative particle production.

Still, the dynamic vacuum can be transported from the baryon or lepton
currents via higher dimensional operators induced from Beyond the SM, which
yields the baryonic or leptonic chemical potential term, like
$\sim\dot{\phi}n_{B/L}$ or $\dot{\phi^{2}}n_{B/L}$, where $n_{B/L}$ denotes
the baryon (lepton) number density and $\phi$ a scalar field which determines
the vacuum. The time-reversal invariance is also locally violated due to the
time-dependent vacuum, so is the C and CP invariance. The original work on the
spontaneous thermal baryogenesis Cohen:1987vi  was based on this possibility,
and recently it has been coupled to the preheating mechanism driven by the
oscillating SM Higgs or the axionlike particles, to be encoded with scenarios
such as the leptogenesis via the axion(like) scalar oscillation
Kusenko:2014uta ; Ibe:2015nfa ; Daido:2015gqa ; Takahashi:2015waa ;
Adshead:2015jza ; Takahashi:2015ula ; Kusenko:2016vcq ; Maleknejad:2016dci ;
DeSimone:2016ofp ; DeSimone:2016juo ; Son:2018avk ; Dasgupta:2018eha ;
Bae:2018mlv ; Domcke:2019qmm ; Wu:2020pej ; Berbig:2023uzs , and the Higgs
relaxation or inflation Kusenko:2014lra ; Yang:2015ida ; Pearce:2015nga ;
Gertov:2016uzs ; Kawasaki:2017ycl ; Wu:2019ohx ; Lee:2020yaj ; Cado:2021bia .
All those scenarios proposed so far have assumed the preheating at higher
scales with the prediction of new particles heavier than the SM particles.

In this paper, we explore a lower-scale preheating with a dynamic and
nonadiabatic baryon chemical potential: it is triggered at the QCD scale of 1
GeV. This is sort of what is called the QCD preheating, which has been
recently proposed in the literature Wang:2022ygk . The presently addressed
scenario is thus alternative to the original one, in a sense that it is the
case when the QCD-induced baryon chemical potential is assumed to be sizable.

The QCD preheating with the dynamical baryon chemical potential is supposed to
be more likely to take place, if the particle production in the universe would
be nonadiabatic. The reason consists of two ingredients: i) as has been
emphasized in Wang:2022ygk , a scalar condensate is present naturally in QCD,
that is the light quark condensate $\sigma\sim\langle\bar{q}q\rangle$, and it
can couple to the nucleon state as well as meson states in a systematic way
respecting the chiral (or isospin) symmetry and its breaking by the current
quark masses; ii) QCD potentially generates the coupling between the dynamic
$\sigma$ field to the baryon current in the chiral (or isospin) symmetric way.
This is indeed the setup that is assumed in higher-scale preheating scenarios
with beyond the SM contents or higher dimensional operators induced from some
beyond the SM. It is all built-in in the QCD preheating case.

An instantaneous out-of-equilibrium can be achieved if QCD is coupled to a
dark QCD (dQCD) sector exhibiting the dark chiral phase transition of the
first order, so that the chiral phase transition in QCD can be supercooled. In
that case the QCD sigma meson field (as the chiral order parameter, or the
light quark condensate) starts to roll in a nonadiabatic way down to the true
QCD vacuum. The induced dynamic baryonic chemical potential is governed by the
dynamic motion of the QCD sigma meson field, so the preheating takes place
along with the dynamic baryon chemical potential: out-of-equilibrium and CP
violation are realized.

We assume QCD to be further allowed to communicate with a dark fermion with
mass of order of 1 GeV and the baryon number violating coupling to neutron,
then the QCD preheating produce the BNA. We discuss this scenario in details
to find that the QCD-induced dynamic baryon chemical potential plays a
significant role for the QCD preheating and the baryogenesis. We show that the
desired amount of the asymmetry today is generated consistently with current
astrophysical, cosmological, and terrestrial experimental constraints.

We also discuss and give comments on the cosmological and phenomenological
consequences characteristic to the present scenario, such as the impacts on
the QCD hadron physics, the cosmological and collider physics probes for the
dQCD sector, and productions of the gravitational waves arising from the
supercooled dQCD phase transition, including comparison between the present
and previous QCD preheating scenarios.

This paper is structured as follows. In Sec. II, we start with introducing a
low-energy effective model of QCD including the QCD sigma meson, pions,
nucleons, and the induced dynamic chemical potential operator, which will be
preliminary setups for the later discussions. In Sec. III we first demonstrate
that thermal QCD baryogenesis, with assuming some Beyond the SM violating the
baryon number symmetry, yields too small BNA. Then, in Sec. IV, we turn to the
preheating scenario driven by the supercooling in the dQCD sector coupled to
QCD. Cosmological and phenomenological constraints on presence of the dQCD
sector are addressed there. The detailed setups relevant to the preheating
analysis are presented in Sec. V. We show the numerical results on the QCD
preheating with the dynamic baryon chemical potential, in Sec. VI, and see
that the desired size of the BNA can be produced due to the QCD-induced
nonadiabatic baryon chemical potential. Section VII provides the summary of
the present paper and several discussions related to the future prospected
probes of the present preheating scenario. Supplemental formulae related to
the numerical analysis on the preheating are listed in Appendices A and B.

##  II The LSM with induced dynamic baryon chemical potential

In this section we start with introducing a low-energy effective description
of QCD including the QCD sigma meson, pions, nucleons, and the induced dynamic
chemical potential operator, which setups the preliminaries in light of the
discussions in the later sections. We employ the linear sigma model (LSM) with
the lightest two flavors $q=(u,d)^{T}$ to monitor the low-energy chiral
dynamics in QCD, based on the chiral $SU(2)_{L}\times SU(2)_{R}$ symmetry and
its breaking structure. The building blocks are: two-by-two complex scalar
matrix: $M\sim\bar{q}_{R}q_{L}$, which transforms under the chiral symmetry as
$M\to U_{L}\cdot M\cdot U_{R}^{\dagger}$ with $U_{L,R}\in SU(2)_{L,R}$, where
$M$ is parametrized by the isosinglet sigma ($\sigma$) mode and isotriplet
pion ($\pi$) mode as $M=\sigma\cdot 1_{\bf 2\times 2}/2+i\pi^{a}\tau^{a}/2$
with the Pauli matrices $\tau^{a}$ ($a=1,2,3$); nucleon (proton,
neutron)-doublet field $N_{L,R}=(p,n)^{T}_{L,R}$, which belong to the
fundamental representation of $SU(2)_{L,R}$ groups.

The LSM Lagrangian is thus given as

$\displaystyle\mathcal{L}_{\rm
LSM}=\operatorname{tr}\left[\partial_{\mu}M^{\dagger}\partial^{\mu}M\right]-V$
$\displaystyle\hskip
20.0pt+\bar{N}i\gamma^{\mu}\partial_{\mu}N-\frac{2m_{N}}{f_{\pi}}\left(\bar{N}_{L}MN_{R}+\bar{N}_{R}M^{\dagger}N_{L}\right)\,,$
$\displaystyle V=m_{\pi}^{2}f_{\pi}\operatorname{tr}\left[{\rm
Re}(M)\right]+m^{2}\operatorname{tr}\left[M^{\dagger}M\right]+\lambda\left(\operatorname{tr}\left[M^{\dagger}M\right]\right)^{2}\,.$
(1)

All the terms are invariant under the chiral symmetry, C and P symmetries,
except for the first term in $V$ coupled to the pion mass, which breaks the
chiral symmetry in a way reflecting the current quark mass term in the
underlying QCD. The explicit-chiral or -isospin breaking for the nucleon
sector has been neglected for simplicity.

The light quark condensate $\langle\bar{q}q\rangle$ generated in QCD is
aligned to be two-flavor diagonal (isospin symmetric) $SU(2)_{V}$ and P
invariant, so that $\langle-\bar{q}_{i}q_{j}\rangle\sim\Lambda_{\rm
QCD}^{3}\cdot\delta_{ij}$, which triggers the (approximate) spontaneous chiral
symmetry breaking $SU(2)_{L}\times SU(2)_{R}\to SU(2)_{V}$ at the vacuum.
Hence the dynamics of the vacuum expectation value of $\sigma$,
$\langle\sigma\rangle$, monitors the dynamic $\langle\bar{q}q\rangle$.

In Eq.(1) and throughout the present study, we take the pion decay constant
$f_{\pi}\simeq 92.4$ MeV, the (isospin symmetric) pion mass $m_{\pi}\simeq
140$ MeV, and the (isospin symmetric nucleon mass) $m_{N}\simeq 940$ MeV. We
fix the potential parameters $m^{2}$ and $\lambda$ in $V$ to the values
determined at the true vacuum $\langle\sigma\rangle=f_{\pi}$ satisfying the
stationary condition $\partial V/\partial\sigma|_{\sigma=f_{\pi}}=0$. Then we
have $\lambda f_{\pi}^{2}=(m_{\sigma}^{2}-m_{\pi}^{2})/2$ and
$m^{2}=\frac{3}{2}m_{\pi}^{2}-\frac{1}{2}m_{\sigma}^{2}$, with the $\sigma$
mass squared defined as
$m_{\sigma}^{2}=\partial^{2}V/\partial\sigma^{2}|_{\sigma=f_{\pi}}$, which we
take
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\simeq(500\>{\rm
MeV})^{2}}$ as the reference mass value of $f_{0}(500)$.

The vacuum expectation value of $\sigma$, $\langle\sigma\rangle$, and its
dynamical evolution are identical to those for $\langle\bar{q}q\rangle$. In
particular, the dynamic $\langle\sigma\rangle$ controls the mass of the
nucleons through the Yukawa interaction in Eq.(1):

$\frac{2m_{N}}{f_{\pi}}\bar{N}_{L}MN_{R}=\frac{m_{N}}{f_{\pi}}\langle\sigma\rangle\cdot\bar{N}_{L}N_{R}+\cdots.$
(2)

Therefore, the nucleon mass
$\tilde{m}_{N}(t)=m_{N}\frac{\langle\sigma(t)\rangle}{f_{\pi}}$ varies in
time, following the time evolution of $\langle\sigma\rangle$.

We further consider the following higher dimensional interaction, which can be
induced from QCD and plays a role of the dynamic baryon chemical potential:

$\displaystyle{\cal L}_{\mu_{\rm dyn}}=-\frac{c}{(4\pi
f_{\pi})^{2}}\cdot\partial_{\mu}{\rm tr}[M^{\dagger}M]J^{\mu}_{B}\,,$ (3)

where $J_{B}^{\mu}$ is the $U(1)$ baryon number current,
$J_{B}^{\mu}=\bar{N}\gamma^{\mu}N$, and the coupling constant $c$ is of ${\cal
O}(1)$ based on the naive dimensional analysis #1#1#1 Since the chemical
potential operator is chiral and QCD gauge invariant, and constructed from the
conserved current, it will not get renormalized, hence no nontrivial
enhancement or suppression of the operator will be generated at the low-energy
scale $\lesssim 4\pi f_{\pi}\sim 1$ GeV. . In the spatially homogeneous
universe, ${\cal L}_{\mu_{\rm dyn}}$ leaves only an effective chemical
potential term dependent of $M(t)$:

$\displaystyle{\cal L}_{\mu_{\rm dyn}}=-\mu_{\rm dyn}(t)\cdot
n_{B}(t)\,,\qquad\mu_{\rm dyn}(t)\equiv\frac{c}{(4\pi
f_{\pi})^{2}}\partial_{0}{\rm
tr}[M(t)^{\dagger}M(t)]=\frac{c}{32\pi^{2}f_{\pi}^{2}}\frac{d}{dt}\langle\sigma^{2}(t)\rangle\,,$
(4)

where $n_{B}(t)=N^{\dagger}(t)N(t)$ is the baryon number density. When the
dynamic $\mu_{\rm dyn}(t)$ gets activated, nonzero dynamic $\mu_{\rm dyn}(t)$
breaks the C and CP symmetries as well as the time-reversal $T$ symmetry. The
form of the dynamic chemical potential operator $\sim\partial_{0}{\rm
tr}[\langle M^{\dagger}M\rangle]$, not including the Nambu-Goldstone boson
mode (i.e. pions), is in contrast to the so-called spontaneous baryogenesis
scenario Cohen:1987vi , or the leptogenesis via the axion(like) scalar
oscillation Kusenko:2014uta ; Ibe:2015nfa ; Daido:2015gqa ; Takahashi:2015waa
; Adshead:2015jza ; Takahashi:2015ula ; Kusenko:2016vcq ; Maleknejad:2016dci ;
DeSimone:2016ofp ; DeSimone:2016juo ; Son:2018avk ; Dasgupta:2018eha ;
Bae:2018mlv ; Domcke:2019qmm ; Wu:2020pej ; Berbig:2023uzs , which is rather
analogous to a class of the leptogenesis via the Higgs relaxation or inflation
Kusenko:2014lra ; Yang:2015ida ; Pearce:2015nga ; Gertov:2016uzs ;
Kawasaki:2017ycl ; Wu:2019ohx ; Lee:2020yaj ; Cado:2021bia .

##  III Thermal QCD-crossover baryogenesis

First of all, we shall consider the ordinary QCD thermal history, in as sense
that $\langle\sigma\rangle\sim\langle\bar{q}q\rangle$ starts to smoothly roll
from $\langle\sigma\rangle\sim 0$ (in the (nearly) chiral symmetric phase)
down to the true vacuum at $\langle\sigma\rangle=f_{\pi}$ during the chiral
crossover transition epoch, $0.7T_{\rm pc}\lesssim T\lesssim T_{\rm pc}$ with
the pseudo critical temperature $\sim 155$ MeV Aoki:2009sc . Once a sufficient
amount of the baryon number violation is supplied from a dark sector beyond
the SM, the thermal QCD baryogenesis would be processed via the dynamic
$\langle\sigma(t)\rangle$ $(0.7T_{\rm pc}\lesssim T\lesssim T_{\rm pc}$) until
$\sigma(t)$ reaches the true vacuum $\langle\sigma\rangle=f_{\pi}$. Here we
assume that the time variation of $\mu_{\rm dyn}(t)$, hence the velocity
$\dot{\sigma}(t)$ is slow enough that the baryon number violating process can
still maintain the thermal equilibrium until $t=t_{D}$, and, when decoupled at
$t=t_{D}$, yield the net asymmetry supplied safely via the thermal chemical
potential $\mu_{\rm dyn}(t_{D})$.

The dynamic chemical potential $\mu_{\rm dyn}(t)$ thus creates the initial net
baryon number $n_{B}(0)$ arising via the thermal distributions for the
nonrelativistic nucleon $N$ and anti-nucleon $\bar{N}$, which can be
estimated, to the first nontrivial order in expansion with respect to
$\mu_{\rm dyn}/T\ll 1$, as

$\displaystyle n_{B}(0)$ $\displaystyle\equiv n_{N}-n_{\bar{N}}$
$\displaystyle=\frac{g_{N}}{2\pi^{2}}\int_{m_{N}}^{\infty}dE\,E(E^{2}-m_{N}^{2})^{1/2}\left[\frac{1}{1+e^{\frac{E-\mu_{\rm
dyn}}{T}}}-\frac{1}{1+e^{\frac{E+\mu_{\rm dyn}}{T}}}\right]$
$\displaystyle\simeq\frac{2g_{N}}{(2\pi)^{3/2}}(m_{N}T)^{3/2}e^{-m_{N}/T}\sinh{\frac{\mu_{\rm
dyn}(t)}{T}}\,,$ (5)

where $g_{N}$ denotes the degrees of freedom for the nucleon. The initial
thermal net number $n_{B}(0)$ in Eq.(5) goes to zero, as $\mu_{\rm dyn}\to 0$.
However, if the baryon number violating process decouples from the thermal
equiliburium faster than $\mu_{\rm dyn}(t)$ reaches zero (i.e.
$\langle\sigma(t)\rangle$ reaches the true vacuum at $f_{\pi}$), the net
number is frozen out at the decoupling temperature $T_{D}$, to be left today.

In the case of the spontaneous baryogenesis, a pseudo Nambu-Goldstone boson
(coupled to the spontaneously broken baryon current) is assumed to be light so
that the mass (i.e., inverse of the damping oscillation time scale) $\ll H$, —
the temperature of the Universe is assumed to be much higher than QCD scale as
well — hence the damping oscillation dynamics of the so-called thermalion
significantly contributes to the thermal baryogenesis (critically assuming the
thermalion to be initially trapped at a slope of the potential) Cohen:1987vi .
In contrast, the thermal QCD phase transition is crossover as observed in
lattice simulations, so that $\langle\sigma(t)\rangle$ is simply expected to
smoothly travel the vacuum transition from $\sim 0$ (at $T=T_{\rm pc}$) to
$f_{\pi}$ (at $T\sim 0.7T_{\rm pc}$) which evolves along with the Hubble
expansion in an adiabatic manner. In that case there is no drastic trapping at
(around) the origin of the potential in contrast to the literature
Cohen:1987vi . The adiabatic motion of $\langle\sigma(t)\rangle$, which
evolves along with the Hubble background, is thus relevant in the present
scenario.

The number density per comoving volume $Y_{B}=n_{B}/s$ is adiabatically
conserved until today: $n_{B}(0)/s(T_{D})=n_{B}(T_{0})/s(T_{0})$ with the
thermal entropy densities $s(T)=\frac{2\pi^{2}}{45}g_{*s}(T)T^{3}$. We refer
to the crossover data on $\langle\sigma(T)\rangle$ from the lattice QCD for
$0.7T_{\rm pc}\lesssim T\lesssim T_{\rm pc}$ Aoki:2009sc  and choosing $T_{D}$
in the crossover regime: $T_{D}={\cal O}(\langle\sigma\rangle)\lesssim{\cal
O}(f_{\pi})$ for $0.7T_{\rm pc}\lesssim T\lesssim T_{\rm pc}$,
$\langle\dot{\sigma}\rangle={\cal O}(f_{\pi}^{2})$, hence $\mu_{\rm dyn}={\cal
O}(c\cdot f_{\pi}/(4\pi)^{2})$ and $\mu_{\rm dyn}/T_{D}={\cal O}(c\cdot
10^{-3})$. Thus we evaluate the yield $Y_{B}(T_{0})$ as

$\displaystyle Y_{B}(T_{0})$ $\displaystyle\sim
10^{-7}\times\left(\frac{\langle\sigma(T_{D}/f_{\pi})\rangle}{f_{\pi}}\right)\cdot\left(\frac{\langle\dot{\sigma}(T_{D}/f_{\pi})\rangle}{f_{\pi}^{2}}\right)\cdot\left(\frac{c}{1}\right)$
$\displaystyle\sim
10^{-7}\times\left(\frac{\langle\sigma(T_{D}/f_{\pi})\rangle}{f_{\pi}}\right)\cdot\left(-\frac{d\langle{\sigma}(T_{D}/f_{\pi})\rangle}{dT}\frac{\frac{T_{D}^{2}}{M_{p}}\cdot
T_{D}}{f_{\pi}^{2}}\right)\cdot\left(\frac{c}{1}\right)$ $\displaystyle\sim
10^{-7}\times\left(\frac{\langle\sigma(T_{D}/f_{\pi})\rangle}{f_{\pi}}\right)\cdot\left(-\frac{d\langle{\sigma}(T_{D}/f_{\pi})\rangle}{dT}\frac{f_{\pi}}{M_{p}}\right)\cdot\left(\frac{c}{1}\right)$
$\displaystyle\sim
10^{-26}\times\left(\frac{\langle\sigma(T_{D}/f_{\pi})\rangle}{f_{\pi}}\right)\cdot\left(\frac{\frac{-d\langle{\sigma}(T_{D}/f_{\pi})\rangle}{dT}}{{\cal
O}(1)}\right)\cdot\left(\frac{c}{1}\right)\,.$ (6)

where , we have roughly read $d[\Delta_{l,s}(T)/\Delta_{l,s}(T=0)]/d[T/T_{\rm
pc}]|_{T\sim T_{\rm pc}\sim T_{D}={\cal O}(f_{\pi})}\sim-{\cal O}(1)$ from the
lattice data on the (normalized) subtracted quark condensate
$\Delta_{l,s}(T)/\Delta_{l,s}(T=0)$ versus $T/T_{\rm pc}$ Aoki:2009sc . This
thermal abundance is too small to be of the order of $10^{-10}$ observed
present day. This is essentially due to the too slow time evolution of
$\langle\dot{\sigma}\rangle$ along with the Planck scale suppression,
following the adiabatic Hubble evolution.

##  IV QCD with dark QCD: triggering the QCD preheating

In the previous section III we have clarified that the QCD-induced dynamic
chemical potential $\mu_{\rm dyn}$ does not work for the thermal baryogenesis,
due to the too slow adiabatic variation of $\mu_{\rm dyn}$, or equivalently,
$\langle\sigma(t)\rangle$ in the chiral crossover transition. In this section
we assume that QCD couples to a dQCD sector with the dark chiral phase
transition of the first order, which makes the QCD chiral phase transition
supercooled and triggers a nonadiabatic roll-down of
$\langle\sigma(t)\rangle$.

We consider $SU(N_{d})$ group for the dark color symmetry of dQCD and
introduce one Dirac fermion $Q_{L,R}$ having ordinary QCD colors ($SU(3)_{c}$)
in the fundamental representation as well: $Q_{L,R}\sim(N_{d},3)$ for
$SU(N_{d})\times SU(3)_{c}$. We set the intrinsic scale of dQCD
$\Lambda_{d}\gtrsim\Lambda_{\rm QCD}\sim 1$ GeV. At the scales below
$\Lambda_{d}\gtrsim\Lambda_{\rm QCD}$, the dark-global chiral
$SU(3)_{dL}\times SU(3)_{dR}$ symmetry will be spontaneously broken by the
dark quark condensate $\langle\bar{Q}Q\rangle$ down to $SU(3)_{dV}$. At almost
the same moment, ordinary QCD will also be confined and generate the
spontaneous breaking of the QCD chiral symmetry via the ordinary quark
condensate $\langle\bar{q}q\rangle$. The low-energy description will thus form
both the ordinary and dark hadron phase reflecting the double confinement as
well.

Here we give a couple of comments on the impact and sensitivity on the hadron
phenomenology in QCD and collider physics probes to be tested in experiments
today. As it turns out, the present dual QCD scenario survives the current
terrestrial experiments.

  * •

At the scale $\Lambda_{d}$, the QCD color nonet
$M_{Q}^{ab}\sim\bar{Q}_{R}^{a}Q_{L}^{b}$ (with $a,b$ being the QCD color
indices) may couple to the QCD quarks and gluons in the gauge- and QCD chiral
invariant way, like

$\displaystyle\bar{q}i\gamma^{\mu}D_{\mu}q\cdot{\rm
tr}[D^{\mu}M_{Q}M_{Q}^{\dagger}]\,,\qquad{\rm
tr}[\bar{q}_{L}q_{R}q_{R}q_{L}]\cdot{\rm
tr}[M_{Q}^{\dagger}M_{Q}]\,,\qquad(\bar{q}i\gamma^{\mu}D_{\mu}q)^{2}\cdot{\rm
tr}[M_{Q}M_{Q}^{\dagger}]\,,\qquad\cdots\,,$ (7)

where ellipses denote terms of much higher dimensional form. We see that the
nonet couplings to quarks do not contribute to the minimal form of the quark-
gluon gauge interaction. Thus the extra hadronic state will necessarily be
exotic to be like eight-quark bound states $\sim qqqqQQQQ$, protected by the
dQCD and QCD color and chiral symmetries. Those exotic hadrons might be
produced via the $\sigma\sigma$ scattering, $\sigma\sigma\to{\rm
tr}[M_{Q}^{\dagger}M_{Q}]\to\sigma\sigma$ including the portal interaction of
the second type in Eq.(7), where the extra resonant signal would be at $\sim
2m_{M_{Q}M_{Q}}\sim 2\Lambda_{\rm QCD}\sim{\rm 2GeV}$ in the invariant mass of
two $\sigma$ mesons.

  * •

Extra light quarks contribute to the running evolution of the ordinary QCD
coupling $g_{s}$, where in the present case $N_{d}$ species of new quarks in
the fundamental representation of $SU(3)_{c}$ group come into play. To keep
the asymptotic freedom, at least $N_{d}$ has to be $<(33/2)-6\sim 10$, which
is determined by the one-loop perturbative calculation. Collider experiments
have confirmed the asymptotic freedom with high accuracy in a wide range of
higher energy scales, in particular above 10 GeV, over 1 TeV
ParticleDataGroup:2022pth . When $\alpha_{s}$ is evolved up to higher scales
using $\alpha_{s}(M_{Z})$, where $\alpha_{s}$ is the fine-structure constant
for $g_{s}$, $\alpha_{s}\equiv g_{s}^{2}/(4\pi)$, measured at the $Z$ boson
pole ($\alpha_{s}(M_{Z})$), as input, (for a recent review of $\alpha_{s}$,
e.g., see Deur:2023dzc ). the tail of the asymptotic freedom around ${\cal
O}(1)$ TeV can thus have sensitivity to exclude new quarks.

Current data on $\alpha_{s}$ at the scale around ${\cal O}(1)$ TeV involve
large theoretical uncertainties. This results in uncertainty of determination
of $\alpha_{s}(M_{Z})$ for various experiments (LHC-ATLAS, -CMS, and Tevatron-
CDF, and D0, etc.), which yields $\alpha_{s}(M_{Z})\simeq 0.110-0.130$, being
consistent with the world average $\alpha_{s}(M_{Z})\simeq 0.118$ within the
uncertainties ParticleDataGroup:2022pth .

We work on the two-loop perturbative computation of $\alpha_{s}$. The dQCD
running coupling ($\alpha_{d}$) contributes to the running of $\alpha_{s}$ at
two-loop level. This contribution is, however, safely negligible: when
$\alpha_{s}\sim\alpha_{d}$ at low-energy scale as desired, because
$\alpha_{d}\ll\alpha_{s}$ at high energy due to the smaller number of dQCD
quarks (with the net number 3 coming in the beta function of $\alpha_{d}$)
than that of the ordinary QCD quarks (with the net number 5 or 6 + $N_{d}$ in
the beta function of $\alpha_{s}$).

Taking into account only the additional $N_{d}$ quark loop contributions to
the running of $\alpha_{s}$, we thus compute the two-loop beta function, and
find that as long as $N_{d}$ is moderately large $(N_{d}\leq 5)$, the measured
ultraviolet scaling (for the renormalization scale of $\mu$ = 10 GeV $-$ a few
TeV) can be consistent with the current data ParticleDataGroup:2022pth  within
the range of $\alpha_{s}(M_{z})$ above. Scaling down to 10 GeV, a prominent
suppression will be seen (e.g. $\alpha_{s}(1\,{\rm GeV})\sim 0.2$ for
$N_{d}=5$ at two-loop level, compared to the perturbative-standard model
prediction $\simeq 0.3$).

Precise measurements in lower scales $\lesssim 10$ GeV have not well been
explored so far, due to the deep infrared complexity of QCD. The low-energy
running of $\alpha_{s}$ is indeed still uncertain, and can be variant as
discussed in a recent review, e.g., Deur:2023dzc . The present dQCD could
dramatically alter the infrared running feature of $\alpha_{s}$, due to new
quarks and the running of the dQCD coupling $\alpha_{d}$. This will also
supply a decisive answer to a possibility of the infrared-near conformality of
real-life QCD.

  * •

Other stringent bounds on the extra light quarks or colored scalars come from
the ALEPH search for gluino and squark pairs tagged with the multijets at the
LEP experiment ALEPH:1997mcm . However, this limit has no sensitivity below
the mass $\sim 2$ GeV, hence is not applicable to the present benchmark model.

Thus a few massless new quarks can still survive constraints on $\alpha_{s}$
at the current status #2#2#2 A similar observation has been made in
Cui:2021sqx  in a different context. . More precise measurements of
$\alpha_{s}$ in the future will clarify how many light or massless new quarks
can be hidden in QCD, which will fix the value of $N_{d}$ in terms of the
present dQCD. The present our aim to consider dQCD is mainly at the dark
chiral phase transition of the first order, which is expected to be realized
in a QCD-like dQCD theory with $N_{d}=3$ and single dark Dirac-fermion flavor
($Q$) having the ordinary QCD colors, which counts the dark flavors to be
three in total. The three-flavor masless QCD-like theory has been supported to
exhbit the first-order chiral phase transition, since the Pisarski-Wilczek’s
original work Pisarski:1983ms  based on the perturbative renormalization group
analysis and also recent functional renormalization group analysis
Fejos:2024bgl . This setup is fully consistent with the current status on the
experimental bound on the extra light quarks as above.

Now we discuss the impact of the doubly QCD charged dQCD quarks on a
supercooled QCD chiral phase transition. In the double confining vacuum, the
QCD-LSM field $M\sim q_{L}q_{R}^{c}$ mixes with the dQCD-QCD color singlet
component $M_{d}\sim\bar{Q}_{L}^{a}Q_{R}^{a}\sim\sigma_{d}+i\eta_{d}$, which
arises from a (nonperturbative) double QCD gluon exchange process in the
ordinary and dark quark scattering with hadronization for the initial and
final states taken into account: in this sense, the process would involve the
mixing with the QCD glueball state $\sim G_{\mu\nu}^{2}$. In the double chiral
invariant manner, it is generated to take the form

$\displaystyle V_{MM_{d}}=\lambda_{\rm mix}{\rm
tr}[M^{\dagger}M](M_{d}^{\dagger}M_{d})\,.$ (8)

Since the generation mechanism is totally nonperturbative, the size of
$\lambda_{\rm mix}$ is expected to be of ${\cal O}(1)$, or so. The dQCD sector
with $N_{d}=3$ in total possesses 3 ($=1_{\rm dark-quarkflavor}\times 3_{\rm
QCD-color}$) flavors, so it is expected to undergo the chiral phase transition
of first order and the critical temperature is $T_{d}={\cal O}(100)$ MeV
$(\sim T_{\rm pc})$. This observation is supported from the Pisarski-Wilczek’s
perturbative renormalization group argument Pisarski:1983ms  and also recent
functional renormalization group analysis Fejos:2024bgl .

Figure 1:  The $\sigma$ potential view schematically drawing different
cosmological setups between the smooth chiral crossover (panel a) and
supercooled chiral phase transition triggered by dQCD (panel b). In the lower
panels a’ and b’ the corresponding thermal evolution of $\sigma$ have also
been displayed. In the panels b and b’, the reheating temperature is evaluated
as $T_{\rm reh}\sim(0.7-0.8)T_{pc}$ as in the text (Eq.(16)), at which epoch
QCD comes back to the thermal equilibrium passing through the out of
equilibrium epoch illustrated by wavy curves in the center of panel b’. More
detailed descriptions are discussed in the text.

Thus at $T=T_{d}$, the chiral order parameter, like $\sigma_{d}^{0}$, drops
down into zero. This is a supercooling, which stores the dark-QCD sector
vacuum energy at the false vacuum. After tunneling from the false vacuum into
the true vacuum, the stored vacuum energy is released into radiations (via the
QCD interaction) and might slightly reheat the universe, keeping the universe
in the thermal equilibrium with $T={\cal O}(T_{\rm pc})$ MeV. In the
literature Sagunski:2023ynd , the supercooled QCD chiral phase transition,
triggered by some dark sector contribution, has been addressed based on the
Nambu-Jona-Lasinio description of low-energy QCD with light three flavors,
where in the case of the vacuum energy as small as the QCD scale, the
percolation temperature has been estimated to be almost identical to the
critical temperature of the first order phase transition. This feature may be
applicable also to the presently assumed dQCD case with three flavors (in
total). Thus the latent heat generated from the dQCD chiral phase transition
would be small enough that the universe can still keep almost the same
temperature of ${\cal O}(T_{\rm pc})$.

Meanwhile, the QCD $\sigma$ gets a discontinuous drop by via the $\lambda_{\rm
mix}$ portal coupling, which is instantaneous and much faster than the Hubble
rate evolving along with the QCD thermal plasma. In the present study, we
assume that the dQCD $\sigma_{d}$ reaches the true vacuum immediately after
the supercooling ends, so that the QCD $\sigma$-fast roll dynamics solely
processes decoupled from the dQCD sector. The $\sigma$ dropping can be
observed in a drastic change of the QCD $\sigma$ mass term at $T=T_{d}$ MeV
due to the induced mass squared $\Delta m^{2}_{\sigma}=\lambda_{\rm
mix}f_{\pi_{d}}^{2}$ with the dark pion decay constant
$f_{\pi_{d}}\sim\langle\sigma_{d}\rangle$, which, when it overwhelms the
thermal mass, drives the discontinuous jump or drop of the QCD $\sigma$, and
triggers the nonadiabatic-fast roll down around $\sigma\sim 0$. Including this
induced $\Delta m^{2}_{\sigma}$ into the original mass square of the QCD
$\sigma$, the net mass parameter realizes the usual QCD vacuum with
$\langle\sigma\rangle=f_{\pi}$. The so-called portal coupling in Eq.(8) thus
plays the important role to trigger the fast roll of $\sigma$ and the QCD
preheating. Figure 1 shows a schematic picture to describe the different
cosmological setups between the chiral crossover and the dQCD-driven
supercooled chiral phase transition.

##  V QCD preheating with nonadiabatic baryon chemical potential

Figure 2:  The $\chi$ lifetime ($\tau_{\chi}$) constraints (black curves) on
the $n-\chi$ mixing coupling $g$ coming from $\tau_{\chi}<\tau_{\rm BBN}\simeq
0.1$s (red straight line), given the mass difference $\Delta m$ ranged from 10
to 200 MeV.

To realize the baryon number violation in the QCD chiral phase transition
epoch, we further consider presence of another dark fermion fields
$\chi_{L,R}$ of the SM singlet, which are allowed to couple to the QCD neutron
fields $n_{L,R}$ in a minimal way as in the literature Bringmann:2018sbs :

$\displaystyle-{\cal L}_{n-\chi}$
$\displaystyle=m_{D}(\bar{\chi}_{R}\chi_{L})+g_{L}(\bar{n}_{R}\chi_{L})+g_{R}(\bar{\chi}_{R}n_{L})+\text{
h.c.}\,.$ (9)

All couplings in Eq.(9) are assumed to be real. The dark Dirac fermion $\chi$
does not carry the baryon number, so $g_{L,R}$ couplings break the baryon
number $(B)$ symmetry as well as the dark fermion number ($B_{\chi}$) symmetry
– but, the total number of $B+B_{\chi}$ is conserved. This type of $n-\chi$
mixing coupling can be generated via leptoquark exchanges involving the
$udd\chi$ vertex, in a way similar to the scenario addressed in the literature
Fornal:2018eol ; Cline:2018ami ; Grinstein:2018ptl ; Keung:2019wpw ;
Keung:2020teb ; Elahi:2020urr ; Fornal:2020bzz ; Fornal:2020gto ;
Fajfer:2020tqf ; Alonso-Alvarez:2021oaj ; Strumia:2021ybk ; Berger:2023ccd ;
Fornal:2023wji ; Liang:2023yta ; Gardner:2023wyl ; Bastero-Gil:2024kjo .

The $\chi$ lifetime is constrained to escape from the BBN and/or CMB
observations: when $m_{\chi}>m_{n}$ as assumed in the present study, $\chi$
having a sizable mixing with the neutron can be completely free from the
existing cosmological and astrophysical constraints, if and only if $\chi$
decays before the BBN epoch at $t_{\rm BBN}\simeq 0.1$s or the lifetime is
longer than the recombination epoch $t_{\rm rec}\sim 10^{13}$s McKeen:2020oyr
. This gives the constraints on the size of the mixing couplings $g_{L,R}$. As
in the case with the literature McKeen:2020oyr , the presently addressed
$\chi$ lifetime is determined also by the transitioned-anomalous magnetic
moment interaction:

$\displaystyle{\cal L}_{\chi
n\gamma}=\frac{\mu_{n}}{2}\cdot\theta\cdot\bar{\chi}\sigma_{\mu\nu}F^{\mu\nu}n+{\rm
h.c.}\,,$ (10)

where $\mu_{n}\simeq 1.91\mu_{N}$ is the neutron magnetic dipole moment with
$\mu_{N}=e/(2m_{n})\simeq 0.1$ e fm being the nuclear magneton
ParticleDataGroup:2022pth . Here we have assumed $g_{L}=g_{R}=g$ for
simplicity, and introduced the mixing angle parameter $\theta=g/\Delta m$ with
the mass difference $\Delta m=m_{\chi}-m_{n}$. The $\chi\to n\gamma$ decay
rate is then computed, for small enough $\theta$, as

$\displaystyle\Gamma[\chi\to n\gamma]\simeq\frac{1}{2200\,{\rm
s}}\left(\frac{\theta}{10^{-10}}\right)^{2}\Bigg{|}\frac{\Delta m}{10\,{\rm
MeV}}\Bigg{|}^{3}\,.$ (11)

Figure 2 shows the lifetime constraint on the mixing coupling $g$ for fixed
$\Delta m$, which implies

$\displaystyle g\gtrsim 10^{-8}\,(10^{-7})\,{\rm MeV}\,,\qquad{\rm
for}\qquad\Delta m=10\,(200)\,{\rm MeV}\,.$ (12)

Combining the mixing coupling term in Eq.(10) with the LSM Lagrangian in
Eq.(1), we work on the following reduced system:

$\displaystyle{\cal L}$
$\displaystyle=\frac{1}{2}(\dot{\sigma})^{2}+\bar{n}\gamma^{\mu}i\partial_{\mu}n+\bar{\chi}\gamma^{\mu}i\partial_{\mu}\chi$
$\displaystyle-\frac{m_{N}}{f_{\pi}}\sigma\bar{n}n-m_{\chi}\bar{\chi}\chi-g(\bar{n}\chi+\bar{\chi}n)+\mu_{\rm
dyn}(t)n^{\dagger}n-V_{\sigma}\,,$ (13)

where

$\displaystyle\mu_{\rm dyn}(t)$
$\displaystyle=\frac{c}{32\pi^{2}}\frac{d}{dt}[{\sigma}^{2}(t)]\,,$
$\displaystyle V_{\sigma}$
$\displaystyle=-m_{\pi}^{2}f_{\pi}\sigma+\frac{1}{2}m^{2}\sigma^{2}+\frac{1}{4}\lambda\sigma^{4}\,,$
(14)

with the same input applied to the potential parameters as in Eq.(1). Here
$\sigma$ field acts as the time-dependent space-homogeneous background field,
$\sigma=\langle\sigma(t)\rangle$, while neutron and dark fermion fields are
quantum fluctuating in spacetime. At this moment the dQCD sector sigma
$\sigma_{d}$ has already reached the true vacuum as noted in the previous Sec.
IV.

The time evolution of the vacuum expectation values in Eq.(37) is evaluated by
solving chained equations of motion for $n$ coupled to $\chi$ through Eq.(9),
together with the dynamic $\langle\sigma(t)\rangle$:

$0=\langle\ddot{\sigma}\rangle+\gamma\langle\dot{\sigma}\rangle-
m_{\pi}^{2}f_{\pi}+m^{2}\langle\sigma\rangle+\lambda\langle\sigma\rangle^{3}+\cdots\,.$
(15)

The ellipses denote the negligible terms including the Hubble friction term
$(3H\langle\dot{\sigma}\rangle)$ and the backreactions from the pion and
nucleon fields #3#3#3 The QCD thermal corrections to the potential of $\sigma$
can also be safely neglected because the potential form gets almost close to
the one at $T=0$, when the roll of $\sigma$ starts at around the reheating
temperature $T\sim 0.7T_{\rm pc}$, as illustrated in Fig. 1. This is supported
from the lattice QCD result Aoki:2009sc  which shows the subtracted light-
quark condensate $\Delta_{l,s}(T)\sim\Delta_{l,s}(T=0)$ at around $T\sim
0.7T_{\rm pc}$. . $\gamma$ plays the role of the full width of the $\sigma$
meson identified as $f_{0}(500)$. As a phenomenological input, we take
$\gamma$ to be the central value of the current measurement,
$\gamma\equiv|2{\rm Im}[\sqrt{s}_{\rm pole}]|=550$ MeV
ParticleDataGroup:2022pth . Here we have applied the conventional vacuum
saturation ansatz: $\langle AB\rangle\equiv\langle A\rangle\langle B\rangle$,
which works fine in the large $N_{c}$ limit. We have also dropped the chemical
potential contribution to the time evolution of $\sigma$ above, because it is
quantitatively small enough due to the loop suppression associated with the
higher-dimensional form of the operator.

Since the dynamic $\langle\sigma\rangle$ instantaneously rolls down from
$\langle\sigma\rangle\sim 0$ to $f_{\pi}$ with the strong friction, the
potential energy
$V(\sigma=0)-V(\sigma=f_{\pi})=\frac{1}{8}(M_{\sigma}^{2}+3m_{\pi}^{2})f_{\pi}^{2}\simeq(135\>{\rm
MeV})^{4}$ is converted into the radiation energy. This triggers reheating of
the universe and the reheating temperature can then be estimated as
Wang:2022ygk

$\displaystyle T_{\rm reh}$
$\displaystyle=\left[\frac{\frac{\pi^{2}}{30}g_{*}(T_{\rm reh})^{4}+(135\>{\rm
MeV})^{4}}{\frac{\pi^{2}}{30}g_{*}}\right]^{1/4}$ $\displaystyle\simeq
119\>{\rm MeV}\sim 0.77T_{pc},$ (16)

where we have assumed the relativistic degrees of freedom $g_{*}(T_{\rm
reh})\simeq 14$ at $T=T_{\rm reh}\sim(0.7-0.8)T_{\rm pc}\sim 110-120\,{\rm
MeV}$ when the dynamic $\langle\sigma\rangle$ starts to roll. The produced
thermal entropy density can also be estimated as

$\displaystyle\frac{2\pi^{2}}{45}g_{*}\cdot(120\>{\rm MeV})^{3}\sim
10^{7}\>{\rm MeV}^{3}.$ (17)

For practical convenience, we reformulate the fermion sector by means of the
two-component spinors. The relations to the Dirac fermions and the definition
of Dirac gamma matrices in terms of the two-component spinor indices are as
follows:

$\displaystyle n$ $\displaystyle=\left(\begin{array}[]{c}n_{L}^{\alpha}\\\
(n_{R}^{c\dagger})^{\dot{\alpha}}\end{array}\right),\qquad\chi=\left(\begin{array}[]{c}\chi_{L}^{\alpha}\\\
(\chi_{R}^{c\dagger})^{\dot{\alpha}}\end{array}\right)\,,$ (22)
$\displaystyle\gamma_{\mu}$
$\displaystyle=\left(\begin{array}[]{cc}0&\sigma_{\mu}^{\alpha\dot{\alpha}}\\\
\bar{\sigma}^{\dot{\alpha}\alpha}&0\end{array}\right)\,,$ (25)

where the superscript “$c$” denotes the charge conjugate. Here
$\alpha,\dot{\alpha}$ are spinor and dotted spinor indices acting on fermion
and anti-fermion fields, and
$\bar{\sigma}_{i}^{\dot{\alpha}\alpha}=-\sigma_{i}^{\dot{\alpha}\alpha}$ for
$i=1,2,3$, $\sigma_{0}^{\alpha\dot{\alpha}}=1^{\alpha\dot{\alpha}}$ and
$\bar{\sigma}_{0}^{\dot{\alpha}\alpha}=1^{\dot{\alpha}\alpha}$. Then Eq.(13)
is cast into the form

$\displaystyle{\cal L}$
$\displaystyle=\frac{1}{2}(\dot{\sigma})^{2}+n^{\dagger}_{L}i\bar{\sigma}^{\mu}\partial_{\mu}n_{L}+n^{c}_{R}i\sigma^{\mu}\partial_{\mu}n_{R}^{c{\dagger}}+\chi^{\dagger}i\bar{\sigma}^{\mu}\partial_{\mu}\chi_{L}+\chi^{c}_{R}i\sigma^{\mu}\partial_{\mu}\chi_{R}^{c{\dagger}}$
$\displaystyle-\frac{m_{N}}{f_{\pi}}\sigma({n}^{\dagger}_{L}n_{R}^{c{\dagger}}+{\rm
h.c.})-m_{\chi}({\chi}^{\dagger}\chi_{R}^{c{\dagger}}+{\rm
h.c.})-g(n_{L}\chi_{R}^{c}+\chi_{L}n_{R}^{c}+{\rm h.c.})$
$\displaystyle+\mu_{\rm
dyn}(t)\left(n_{L}^{\dagger}\bar{\sigma}_{0}n_{L}+n_{R}^{c}{\sigma}_{0}n_{R}^{c{\dagger}}\right)$
$\displaystyle-V\,.$ (26)

The relevant equations of motion for the fluctuating spinor fields are:

$\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}n_{L}+\tilde{m}_{N}n_{R}^{c\dagger}+g\chi_{R}^{c\dagger}+\mu_{\rm
dyn}(t)\bar{\sigma}_{0}n_{L}\,,$ (27) $\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}n_{R}^{c{\dagger}}+\tilde{m}_{N}n_{L}+g\chi_{L}+\mu_{\rm
dyn}(t)\sigma_{0}n_{R}^{c{\dagger}}\,,$ $\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}\chi_{L}+m_{\chi}\chi_{R}^{c\dagger}+gn_{R}^{c\dagger}\,,$
$\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}\chi_{R}^{c{\dagger}}+m_{\chi}\chi_{L}+gn_{L}\,,$

where

$\displaystyle\tilde{m}_{N}(t)=m_{N}\frac{\langle\sigma(t)\rangle}{f_{\pi}}\,.$
(28)

The time-varying mass $\tilde{m}_{N}$ in Eq.(28) causes the nonperturbative
nucleon production when the adiabaticity is violated:
$|\dot{\tilde{m}}_{N}/\tilde{m}_{N}^{2}|\gtrsim 1$. This inequality leads to
the production range in terms of the $\sigma$ motion as

$|\sigma|\lesssim\sqrt{\frac{f_{\pi}\langle\dot{\sigma}\rangle}{m_{N}}}\simeq
42\>{\rm MeV},$ (29)

where we read
$\langle\dot{\sigma}\rangle\simeq\sqrt{V(\sigma=0)-V(\sigma=f_{\pi})}\sim(135\>{\rm
MeV})^{2}$. The actual production time can be earlier because the estimated
velocity $\langle\dot{\sigma}\rangle$ would be smaller due to the friction
$\gamma$.

We move on to the momentum space by the Fourier transform:

$\left(\begin{array}[]{c}n_{L}(t,\vec{x})\\\ n_{R}^{cT}(t,\vec{x})\\\
\chi_{L}(t,\vec{x})\\\
\chi_{R}^{cT}(t,\vec{x})\end{array}\right)_{\alpha}=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}}\sum_{s=\pm}\left(e_{\vec{k}}^{s}\right)_{\alpha}\left(\begin{array}[]{c}n_{L,\vec{k}}^{s}(t)\\\
n_{R,\vec{k}}^{s}(t)\\\ D_{L,\vec{k}}^{s}(t)\\\
D_{R,\vec{k}}^{s}(t)\end{array}\right)\,,$ (30)

where we have introduced the helicity eigenstate vector having helicity
$s=+,-$ as

$\left(e_{\vec{k}}^{s}\right)_{\alpha}=\left(\begin{array}[]{c}\sqrt{\frac{1}{2}\left(1+\frac{sk^{3}}{k}\right)}\\\
s\,e^{i\theta_{k}}\sqrt{\frac{1}{2}\left(1-\frac{sk^{3}}{k}\right)}\end{array}\right)\,,$
(31)

with $k\equiv|\vec{k}|$ and

$e^{i\theta_{\vec{k}}}\equiv\frac{k^{1}+ik^{2}}{\sqrt{\left(k^{1}\right)^{2}+\left(k^{2}\right)^{2}}}\,,$
(32)

which are the solutions to the helicity-eigen equations for left- and right-
handed spinors:

$\displaystyle-\frac{k^{i}}{k}(\sigma_{0}\bar{\sigma}_{i})^{\beta}_{\alpha}\cdot(e_{\vec{k}}^{s})_{\beta}$
$\displaystyle=s\cdot(e^{s}_{\vec{k}})_{\alpha}\,,$
$\displaystyle-\frac{k^{i}}{k}(\bar{\sigma}_{0}\sigma_{i})_{\dot{\beta}}^{\dot{\alpha}}\cdot(e_{\vec{k}}^{s{\dagger}})^{\dot{\beta}}$
$\displaystyle=s\cdot(e^{s{\dagger}}_{\vec{k}})^{\dot{\alpha}}\,.$ (33)

Applying those transformed fields into equations of motion in Eq.(27), we get

$\displaystyle{\cal D}_{t}n_{L,\vec{k}}^{s}$
$\displaystyle=iskn_{L,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}n_{R,-\vec{k}}^{s\dagger}+gD_{R,-\vec{k}}^{s\dagger}\right)\,,$
(34) $\displaystyle{\cal D}_{t}n_{R,\vec{k}}^{s}$
$\displaystyle=iskn_{R,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}n_{L,-\vec{k}}^{s\dagger}+gD_{L,-\vec{k}}^{s\dagger}\right)\,,$
$\displaystyle{\cal D}_{t}D_{L,\vec{k}}^{s}$
$\displaystyle=iskD_{L,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}D_{R,-\vec{k}}^{s\dagger}+gn_{R,-\vec{k}}^{s\dagger}\right)\,,$
$\displaystyle{\cal D}_{t}D_{R,\vec{k}}^{s}$
$\displaystyle=iskD_{R,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}D_{L,-\vec{k}}^{s\dagger}+gn_{L,-\vec{k}}^{s\dagger}\right)\,,$

where

$\displaystyle{\cal D}_{t}n_{L,\vec{k}}^{s}$
$\displaystyle=(\partial_{t}-i\mu_{\rm dyn}(t))n_{L,\vec{k}}^{s}\,,$
$\displaystyle{\cal D}_{t}n_{R,\vec{k}}^{s}$
$\displaystyle=(\partial_{t}+i\mu_{\rm dyn}(t))n_{R,\vec{k}}^{s}\,,$
$\displaystyle{\cal D}_{t}D_{L,\vec{k}}^{s}$
$\displaystyle=\partial_{t}D_{L,\vec{k}}^{s}\,,$ $\displaystyle{\cal
D}_{t}D_{R,\vec{k}}^{s}$ $\displaystyle=\partial_{t}D_{R,\vec{k}}^{s}\,,$ (35)

and we have used the orthogonality and completeness relations for helicity
eigenstates,

$\displaystyle
e^{s{\dagger}}_{\vec{k}}\bar{\sigma}_{0}e^{s^{\prime}}_{\vec{k}}=e^{s}_{\vec{k}}{\sigma}_{0}e^{s^{\prime}{\dagger}}_{\vec{k}}=\delta^{ss^{\prime}}\,,$
$\displaystyle\sum{s,s^{\prime}}(e^{s}_{\vec{k}})_{\alpha}(e^{s{\dagger}}_{\vec{k}})_{\dot{\alpha}}=(\sigma_{0})_{\alpha\dot{\alpha}}\,,\qquad\sum{s,s^{\prime}}(e^{s{\dagger}}_{\vec{k}})^{\dot{\alpha}}(e^{s}_{\vec{k}})^{{\alpha}}=(\bar{\sigma}_{0})_{\dot{\alpha}\alpha}\,,$
$\displaystyle
e^{s}_{\vec{k}}e^{s^{\prime}}_{-\vec{k}}=s\,e^{i\theta_{\vec{k}}}\cdot\delta^{ss^{\prime}}\,,\qquad
e^{s{\dagger}}_{-\vec{k}}e^{s^{\prime}{\dagger}}_{\vec{k}}=s\,e^{-i\theta_{\vec{k}}}\cdot\delta^{ss^{\prime}}\,.$
(36)

Because of the spacial-homogeneous universe we are allowed to focus only on
one mode ($n_{k}$) having $\vec{k}=\left(\epsilon,0,k_{3}\right)$ with
$\epsilon\rightarrow 0$ in the end, which leads to $e^{i\theta_{\vec{k}}}=1$,
and then perform the momentum-space integral $\int 4\pi k^{2}n_{k}dk$.

The net baryon number density is computed as

$\displaystyle n_{B}(t)$ $\displaystyle=\frac{1}{V}\int d^{3}x\left(\langle
n^{\dagger}n\rangle-\sum_{\vec{k}}2\right)\,,$ (37)

where $V=\int d^{3}x$ is the space volume, and the last term corresponds to
the subtraction of the divergent part induced by the zero-point energy
($4\times\frac{1}{2}$). The net baryon number density in Eq.(37) can be
rewritten in terms of the vacuum expectation values of two Fourier transformed
fields as

$\begin{aligned}
n_{B}(t)&=\frac{1}{V}\int\frac{d^{3}k}{(2\pi)^{3}}\sum_{s=\pm}\left(\left\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle-\left\langle
n_{R,\vec{k}}^{s\dagger}n_{R,\vec{k}}^{s}\right\rangle\right)\end{aligned}\,,$
(38)

The two-point functions relevant to this net baryon number density are
developed in time through the time evolution equations given in Eq.(34). The
explicit expressions of those time evolution equations and the initial
conditions for the two-point functions are supplemented in Appendix B.

##  VI Baryogenesis via QCD preheating

In this section, we numerically evaluate the time-evolution of
$\langle\sigma(t)\rangle$, the chemical potential $\mu_{dyn}(t)$, and the BNA
$Y_{B}=\frac{n_{B}}{s}$, based on the formulae presented in the previous Sec.
V and listed in Appendix B. The entropy density $s$ is set as the thermal one,
in Eq.(17), which turns out to be comparable with the nonperturbative pionic
entropy density created through the preheating as noted in Wang:2022ygk .

$\left\langle\sigma(t)\right\rangle$ is evaluated by solely solving the space-
homogeneous (approximate) time-evolution Eq.(15) with the input parameters
setup as in the previous sections. First of all, we recall that the preheating
scalar $\sigma$ starts to roll down from $\sigma\sim 0$, which is driven via
the portal coupling in Eq.(8) instantaneously just after the dQCD $\sigma_{d}$
undergoes the first-order phase transition at $T\sim f_{\pi_{d}}={\cal
O}(100)$ MeV, as in the right panel of Fig. 1. Thus we can simply set the
initial condition of $\langle\sigma(t)\rangle$ as

$\displaystyle\langle\sigma(t=t_{0})\rangle=\langle\dot{\sigma}(t=t_{0})\rangle=0\,.$
(39)

The nonzero $\langle\sigma(t=0)\rangle$ will not substantially alter the
essential features of the nonadiabatic evolution of $\langle\sigma(t)\rangle$.
Recall that the nonadiabatic condition is met when
$|\sigma|\lesssim\sqrt{\frac{f_{\pi}\langle\dot{\sigma}\rangle}{m_{N}}}\simeq
42\>{\rm MeV}$ (Eq.(29)), which corresponds to

$M_{\sigma}\cdot t\lesssim 5\,.$ (40)

See the orange curve in Fig. 3. Since the time scale of the $\sigma$ motion
($\sim M_{\sigma}^{-1}$) is thus much shorter than the Hubble rate $(\sim
M_{p}/T^{2}_{\rm pc})$, the produced neutrons and pions will never be
thermalized.

The total nucleon number density can be evaluated as

$\displaystyle n_{N}(t)+\bar{n}_{N}(t)$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}}\frac{\langle\tilde{\rho}_{N}(\vec{k},t)\rangle}{\omega(\vec{k},t)}\simeq
2\sum_{\vec{k}}\frac{\langle\tilde{\rho}_{n}(\vec{k},t)\rangle}{\omega(\vec{k},t)}\,.$
(41)

for $N=p,n$, where $\vec{A}$ denotes the spatial component of the vector
variable $A$ and $\omega(\vec{k},t)=\sqrt{|\vec{k}|^{2}+\tilde{m}_{N}^{2}(t)}$
is the one-particle energy of the nucleon. In reaching the second equality of
Eq.(41) we have simply taken the isospin symmetric limit.
$\tilde{\rho}_{n}(\vec{k},t)$ is the kinetic energy density of the nucleon in
momentum space, which is derived from the Hamiltonian as #4#4#4 A similar
formula for the energy density in Wang:2022ygk  includes errors, which have
fully been corrected in Eq.(41).

$\displaystyle\langle\tilde{\rho}_{n}(\vec{k},t)\rangle$
$\displaystyle=\frac{1}{V}\sum_{s=\pm}\left[\omega(\vec{k},t)-s|\vec{k}|\left(\left\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle+\left\langle
n_{R,\vec{k}}^{s\dagger}n_{R,\vec{k}}^{s}\right\rangle\right)\right.$
$\displaystyle\left.+s\,\tilde{m}_{N}\left(e^{i\theta_{\vec{k}}}\left\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\right\rangle+e^{-i\theta_{\vec{k}}}\left\langle
n_{L,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\right\rangle\right)\right]\,.$
(42)

The detailed derivation of Eq.(42) is supplemented in Appendix A. The last
term in Eq.(42) corresponds to the subtraction of the negative vacuum energy
($4\times\frac{1}{2}\omega(\vec{k},t)$). The time evolution of the total
baryon number density in Eq.(41) is thus determined by solving coupled
equations of motion for $n_{L,\vec{k}}^{s}$, $n_{R,\vec{k}}^{s}$, and
$\langle\sigma(t)\rangle$ (Eqs.(34) and (15)). Note that as long as the
$n-\chi$ mixing coupling $g$ is small enough compared to the typical scale of
the preheating frequency $\sim{\cal O}(M_{\sigma})$, the total number density
$(n_{N}+\bar{n}_{N})$ is almost completely determined by QCD without dQCD, as
was stressed in Wang:2022ygk . The dynamical chemical potential term has been
also dropped in there, because the total number density $(n_{N}+\bar{n}_{N})$
is dominantly produced via the strong Yukawa coupling $\sim
m_{N}/f_{\pi}={\cal O}(10)$.

$(n_{N}+\bar{n}_{N})$ is explosively generated by the nonadiabatic-time
varying $\langle\sigma(t)\rangle$ and gets asymptotically saturated to be
$\simeq 10^{5}\ {\rm MeV}^{3}$ in the time range $M_{\sigma}t\lesssim 5$
(Eq.(40)). This number density is much larger than the thermal equilibrium
density at the reheating temperature $T_{\rm reh}\simeq 119$ MeV (Eq.(16)),
$[(n_{N}+\bar{n}_{N})]_{\rm EQ}=2\times 2\times 2\times(m_{N}T_{\rm
reh}/2\pi)^{3/2}e^{-m_{N}/T_{\rm reh}}\sim 7\times 10^{3}\ {\rm MeV}^{3}$,
hence becomes overabundant. Actually, the overproduced nucleons can survive
long enough during the whole reheating process: the relaxation time scale, at
which the overproduced nucleons pair-annihilate, can be estimated as

$\displaystyle M_{\sigma}\Delta t_{\rm relax}=M_{\sigma}(n_{N}\langle\sigma
v\rangle)^{-1}\sim 700\,,$ (43)

or equivalently, $\Delta t_{\rm relax}\sim 5\times 10^{-7}\ {\rm fs}$ for
$M_{\sigma=f_{0}(500)}=500$ MeV. Here the static nucleon-pair annihilation
cross section has been evaluated as a classical disc $\langle\sigma
v\rangle\sim 4\pi/m_{N}^{2}$ with the impact parameter $(1/m_{N})$. Thus the
QCD preheating is operative to accumulate a large number of nucleon and
antinucleon pairs in out-of-equilibrium until the relaxation time. This is a
novel picture of the thermal history in the QCD phase transition epoch, as
emphasized in the literature Wang:2022ygk .

Figure 3: The time evolution of $\langle\sigma(t)\rangle$ (orange curve) and
$\mu_{\rm dyn}$ in Eq.(4) (blue curve) in the QCD preheating epoch defined in
the time interval $0\leq M_{\sigma}t\leq 700$ as in Eq.(43).  Figure 4: The
time evolution and freeze out of BNA in the QCD preheating epoch, as in Fig.
3, for various sizes of the $n-\chi$ mixing coupling strength $g$. The dashed
vertical line at $M_{\sigma}t=700$ corresponds to the relaxation time scale in
Eq.(43) until which scale the QCD preheating needs to be finalized. In this
figure we have taken $m_{N}=1000$ MeV, $m_{\chi}=1050$ MeV, and $c=1$, as a
reference point. For the parameter dependence, see the text.

In Fig 3 we show the time evolution of $\langle\sigma(t)\rangle$ (orange
curve) and $\mu_{\rm dyn}(t)$ (blue curve) until around the relaxation time
scale in Eq.(43). $\langle\sigma(t)\rangle$ reaches the true vacuum
$\langle\sigma\rangle=f_{\pi}$ at $M_{\sigma}t\sim 10$ passing through the
nonadiabatic particle production regime $M_{\sigma}t<5$ in Eq.(40). It then
goes through the true vacuum $\langle\sigma\rangle=f_{\pi}$ to be pulled back
and trapped at the true vacuum due to the strong friction of the $\sigma$
decay width $\gamma$. Thus the $\langle\sigma(t)\rangle$ evolution does not
simply go like $\sim\cos M_{\sigma}t$, hence the particle production does not
follow the Mathieu equation or the parametric resonance Wang:2022ygk .

The dynamic baryon chemical potential $\mu_{\rm dyn}$, which essentially
scales with $\langle\dot{\sigma}\rangle$ as in Eq.(4), grows up until
$\langle\sigma(t)\rangle$ reaches the true vacuum, passing through the
nonadiabatic production epoch. It then drastically damps due to $\gamma$, and
finally freezes out at $M_{\sigma}t={\cal O}(100)$ (which keeps exhibiting a
tiny intrinsic oscillation in magnitude, though).

Figure 4 displays the BNA evolved until the relaxation time scale
$M_{\sigma}t\sim 700$ with varying the $n-\chi$ mixing strength $g$ (in unit
of MeV). The BNA yield $Y_{B}$ grows in magnitude until
$\langle\sigma(t)\rangle$ reaches the true vacuum for the first time (at
$M_{\sigma}t\sim 10$ as seen in Fig. 3). The BNA yield $Y_{B}$ is essentially
produced by the nonadiabatic $\sigma$ motion and the associated
nonperturbative nucleon-antinucleon pair production for the time scale in the
production era $0\leq M_{\sigma}t\lesssim 5$. Finally, $Y_{B}$ asymptotically
gets saturated to constant due to the decoupling of $\mu_{\rm dyn}$, which
starts to operate at around $M_{\sigma}t\sim 10$ (see Fig. 3) and drives the
$n-\bar{n}$ oscillation to end until the relaxation time scale
$M_{\sigma}t\sim 700$.

The size of the net baryon number density $n_{B}$ evolves like
$|n_{B}(t)|\sim|\frac{\mu_{\rm dyn}(t)g^{4}}{m_{\chi}^{2}-m_{N}^{2}}|$, so
simply grows with the $n-\chi$ mixing coupling strength $g$ and the size of
the dynamic baryon chemical potential, $c$ in Eq.(4). We have also checked
that the dependence of the dark fermion mass $m_{\chi}$ is small enough that
the BNA time evolution until the relaxation time scale is not substantially
altered from the case with $m_{\chi}=1050$ MeV, as in Fig. 4, as long as
$m_{\chi}$ is close, but greater than $m_{N}$ (such as in a range of
$1000<m_{\chi}<1200$ MeV).

The baryogenesis is completed before the relaxation time period
($M_{\sigma}t\sim 700$) and, with $g={\cal O}(10^{-2})$ MeV, successfully
yield the desired amount of the baryon asymmetry, $Y_{B}=n_{B}/s\sim
10^{-10}-10^{-9}$, consistently with the current cosmological and
astrophysical constraint in Eq.(12). Thus, the nonadiabatic $\mu_{\rm dyn}(t)$
significantly works for the baryogenesis via the QCD preheating.

##  VII Summary and Discussions

In summary, the dynamic chemical potential induced from QCD can provide the
key role to generate the BNA in a nonadiabatic way through the QCD preheating
played by the light quark condensate. The idea of the generation mechanism
with the induced chemical potential is analogous to the original spontaneous
baryogenesis Cohen:1987vi , which was based on the adiabatic (thermal)
generation, and the leptogenesis via the Higgs relaxation or inflation
Kusenko:2014lra ; Yang:2015ida ; Pearce:2015nga ; Gertov:2016uzs ;
Kawasaki:2017ycl ; Wu:2019ohx ; Lee:2020yaj ; Cado:2021bia , the axion(like)
scalar oscillation Kusenko:2014uta ; Ibe:2015nfa ; Daido:2015gqa ;
Takahashi:2015waa ; Adshead:2015jza ; Takahashi:2015ula ; Kusenko:2016vcq ;
Maleknejad:2016dci ; DeSimone:2016ofp ; DeSimone:2016juo ; Son:2018avk ;
Dasgupta:2018eha ; Bae:2018mlv ; Domcke:2019qmm ; Wu:2020pej ; Berbig:2023uzs
. In contrast to the latter cases, the presently addressed scenario does not
assume anything at higher scales beyond the SM for generating the chemical
potential operator, which arises from QCD of the SM, instead, predicts new
physics sector (what we call dQCD and a dark Dirac fermion $\chi$) exhibiting
the cosmological phase transition with the supercooling at around the QCD
scale.

The presently proposed baryogenesis includes two dark sectors: a dark fermion
$\chi$ and dQCD with QCD colors. In closing, we give several comments on the
phenomenological and cosmological probes related to those dark sectors, in
addition to the impacts on the QCD hadron physics and TeV scale probes at the
LHC as noted in Sec.IV.

The dark fermion $\chi$ coupled to the neutron, which violates the baryon
number symmetry, can be present with the mass of $\sim 1$ GeV consistently
with the current cosmological and astrophysical bound (Eq.(12)) and the
desired BNA yield (Fig. 4). This baryon number violating coupling may imply
another new physics at higher scales, like leptoquarks as addressed in the
literature Fornal:2018eol ; Cline:2018ami ; Grinstein:2018ptl ; Keung:2019wpw
; Keung:2020teb ; Elahi:2020urr ; Fornal:2020bzz ; Fornal:2020gto ;
Fajfer:2020tqf ; Alonso-Alvarez:2021oaj ; Strumia:2021ybk ; Berger:2023ccd ;
Fornal:2023wji ; Liang:2023yta ; Gardner:2023wyl ; Bastero-Gil:2024kjo , which
would be constrained by the current experimental limits on searching for
leptoquarks (for a review, e.g., see Dorsner:2016wpm ). Thus the allowed
baryon number violating coupling strength $g$ could be sensitive to those
leptoquark searches, depending on the modeling, which might be challenging to
generate the desired BNA as in Fig. 4 and to even survive the BBN limit in
Eq.(12).

The leptoquark physics would also discriminate the present scenario and the
original QCD preheating setup Wang:2022ygk , in which $\mu_{\rm dyn}(t)$ was
simply assumed to be negligible, instead, both the CP violation and the source
to freeze out the BNA were provided from another dark sector. In particular,
in the original work Wang:2022ygk  an axionlike particle having a time-damping
evolution has been assumed to make the baryon number violating coupling
dropped to zero, so as to end the $n-\bar{n}$ oscillation. Searching for
leptoquarks as well as the axionlike particle could thus be a punchline to
test whether the QCD-induced dynamic baryonic chemical potential can be
sizable and baryon number violating coupling is still present today.

The supercooled dQCD in the present scenario can generate the gravitational
wave spectra associated with the percolation/nucleation temperature around the
QCD scale. The signals could be to be probed via the pulsar timing arrays
sensitive to the peak frequencies around nano hertz NANOGrav:2023gor ;
NANOGrav:2023hfp ; EPTA:2023sfo ; EPTA:2023akd ; EPTA:2023fyk ;
Reardon:2023gzh ; Reardon:2023zen ; Xu:2023wog , or the future prospected
interferometers sensitive to higher frequencies LISA:2017pwj ; Caprini:2019egz
; Corbin:2005ny ; Harry:2006fi ; Kawamura:2006up ; Yagi:2011wg . Consistency
with the realization of the QCD preheating setup, as illustrated in Fig. 1,
can thus be discussed in light of the gravitation wave prediction.

Dark matter candidates can arise from the dQCD sector, such as the dark
baryons and/or mesons. The sizable enough relic abundance could be generated
through another preheating which can arise due to the nonadiabatic oscillation
of the dark sigma meson $\sigma_{d}$, triggered by the dQCD supercooling,
right before the QCD preheating epoch. Since the $\sigma-\sigma_{d}$ portal
coupling is present and the dQCD as well as QCD interactions are strong, the
dark-matter direct detection experiments could probe those dQCD dark-matter
hadrons.

Thus, the presently proposed and explored baryogenesis can involve plenty of
phenomenological and cosmological probes to be tested by the current and
future prospected observations in the sky and/or terrestrial experiments.
Those would be deserved to purse in the future publications.

## Acknowledgments

We are grateful to Seishi Enomoto, Yuepeng Guan, Hiroyuki Ishida, and Jinyang
Li for fruitful discussions. This work was supported in part by the National
Science Foundation of China (NSFC) under Grant No.11747308, 11975108,
12047569, and the Seeds Funding of Jilin University (S.M.).

##  Appendix A Derivation of Eq.(42)

###  A.1 Derivation of Eq.(42)

We consider the relevant part of the Lagrangian in Eq.(26):

${\cal
L}=n_{L}^{\dagger}i\bar{\sigma}^{\mu}\partial_{\mu}n_{L}+n^{c}_{R}i\sigma^{\mu}\partial_{\mu}n_{R}^{c{\dagger}}-\tilde{m}_{N}({n}_{L}^{\dagger}n_{R}^{c{\dagger}}+{\rm
h.c.})\,.$ (44)

From this we read off the canonical conjugate fields,

$\tilde{n}_{L}=\frac{\partial{\cal
L}}{\partial\dot{n}_{L}}=n^{{\dagger}}_{L}i\bar{\sigma}_{0};\ \ \ \
\tilde{n}^{c{\dagger}}_{R}=\frac{\partial{\cal
L}}{\partial\dot{n}^{c{\dagger}}_{R}}=n^{c}_{R}i\sigma_{0}\,.$ (45)

By performing the Legendre transformation, we get the Hamiltonian

$\displaystyle H=$
$\displaystyle\int{d^{3}x}(\tilde{n}_{L}\dot{n}_{L}+\tilde{n}^{c{\dagger}}_{R}\dot{n}^{c{\dagger}}_{R}-{\cal
L})$ (46) $\displaystyle=$
$\displaystyle\int{d^{3}x}[-n_{L}^{{\dagger}}i\bar{\sigma}^{i}\partial_{i}n_{L}-n_{R}^{c}i\sigma^{i}\partial_{i}n_{R}^{c{\dagger}}+\tilde{m}_{N}(n_{L}^{{\dagger}}n_{R}^{c{\dagger}}+h.c.)]\,.$

This can further be evaluated in terms of the Fourier transformed fields in
Eq.(30), as follows:

$\displaystyle H=\int{d^{3}x}$
$\displaystyle[-\int{\frac{d^{3}k}{(2\pi)^{3}}}e^{-i\vec{k}\cdot\vec{x}}\sum_{s}e_{\vec{k}}^{s{\dagger}}n_{L,\vec{k}}^{s{\dagger}}\cdot
i\bar{\sigma}^{i}\int{\frac{d^{3}k^{\prime}}{(2\pi)^{3}}}(i\vec{k}^{\prime})e^{-i\vec{k}^{\prime}\cdot\vec{x}}\sum_{s^{\prime}}e_{\vec{k}^{\prime}}^{s^{\prime}}n_{L,\vec{k}^{\prime}}^{s^{\prime}}$
(47)
$\displaystyle-\int{\frac{d^{3}k}{(2\pi)^{3}}}e^{i\vec{k}\cdot\vec{x}}\sum_{s}e_{\vec{k}}^{s}n_{R,\vec{k}}^{s}\cdot
i\sigma^{i}\int{\frac{d^{3}k^{\prime}}{(2\pi)^{3}}}(-i\vec{k}^{\prime})e^{-i\vec{k}^{\prime}\cdot\vec{x}}\sum_{s^{\prime}}e_{\vec{k}^{\prime}}^{s^{\prime}{\dagger}}n_{R,\vec{k}^{\prime}}^{s^{\prime}{\dagger}}$
$\displaystyle+\tilde{m}_{N}(\int{\frac{d^{3}k}{(2\pi)^{3}}}e^{-i\vec{k}\cdot\vec{x}}\sum_{s}e_{\vec{k}}^{s{\dagger}}n_{L,\vec{k}}^{s{\dagger}}\int{\frac{d^{3}k^{\prime}}{(2\pi)^{3}}}e^{-i\vec{k}^{\prime}\cdot\vec{x}}\sum_{s^{\prime}}e_{\vec{k}^{\prime}}^{s^{\prime}{\dagger}}n_{R,\vec{k}^{\prime}}^{s^{\prime}{\dagger}}+h.c.)]$
$\displaystyle=\int{\frac{d^{3}k}{(2\pi)^{3}}\frac{d^{3}k^{\prime}}{(2\pi)^{3}}}$
$\displaystyle[\delta(\vec{k}-\vec{k}^{\prime})\cdot
k_{i}\sum_{s,s^{\prime}}(e_{\vec{k}}^{s{\dagger}}e_{\vec{k}^{\prime}}^{s^{\prime}}n_{L,\vec{k}}^{s{\dagger}}\bar{\sigma}^{i}n_{L,\vec{k}^{\prime}}^{s^{\prime}}-e_{\vec{k}}^{s}e_{\vec{k}^{\prime}}^{s^{\prime}{\dagger}}n_{R,\vec{k}}^{s}\sigma^{i}n_{R,\vec{k}^{\prime}}^{s^{\prime}{\dagger}})$
$\displaystyle+\tilde{m}_{N}\delta(-\vec{k}-\vec{k}^{\prime})\sum_{s,s^{\prime}}(e_{\vec{k}}^{s{\dagger}}e_{\vec{k}^{\prime}}^{s^{\prime}{\dagger}}n_{L,\vec{k}}^{s{\dagger}}n_{R,\vec{k}^{\prime}}^{s^{\prime}{\dagger}}+h.c.)]$
$\displaystyle=\int{\frac{d^{3}k}{(2\pi)^{3}}}\sum_{s}$
$\displaystyle[k_{i}(n_{L,\vec{k}}^{s{\dagger}}\bar{\sigma}^{i}n_{L,\vec{k}}^{s}-n_{R,\vec{k}}^{s}\sigma^{i}n_{R,\vec{k}}^{s{\dagger}})+\tilde{m}_{N}(se^{i\theta_{\vec{k}}}n_{L,\vec{k}}^{s{\dagger}}n_{R,-\vec{k}}^{s{\dagger}}+se^{-i\theta_{\vec{k}}}n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s})]$
$\displaystyle=\int{\frac{d^{3}k}{(2\pi)^{3}}}\sum_{s}$
$\displaystyle[-s\lvert\vec{k}\rvert(n_{L,\vec{k}}^{s{\dagger}}n_{L,\vec{k}}^{s}+n_{R,\vec{k}}^{s}n_{R,\vec{k}}^{s{\dagger}})+\tilde{m}_{N}(se^{i\theta_{\vec{k}}}n_{L,\vec{k}}^{s{\dagger}}n_{R,-\vec{k}}^{s{\dagger}}+se^{-i\theta_{\vec{k}}}n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s})]\,.$

Here we have used the orthogonality relations:

$\displaystyle e_{\vec{k}}^{s{\dagger}}e_{\vec{k}}^{s^{\prime}}$
$\displaystyle=\delta^{ss^{\prime}}\,,$ $\displaystyle
e_{\vec{k}}^{s}e_{\vec{k}}^{s^{\prime}{\dagger}}$
$\displaystyle=\delta^{ss^{\prime}}\,,$ $\displaystyle
e_{\vec{k}}^{s{\dagger}}e_{-\vec{k}}^{s^{\prime}{\dagger}}$
$\displaystyle=\delta_{ss^{\prime}}\cdot se^{i\theta_{\vec{k}}}\,,$
$\displaystyle e_{-\vec{k}}^{s}e_{\vec{k}}^{s^{\prime}}$
$\displaystyle=\delta_{ss^{\prime}}\cdot se^{-i\theta_{\vec{k}}}\,.$ (48)

Thus the kinetic energy density operator, normalized by the volume factor $V$,
is read off as

$\rho_{\vec{k}}=\frac{1}{V}\sum_{s=\pm}[-s\lvert\vec{k}\rvert(n_{L,\vec{k}}^{s{\dagger}}n_{L,\vec{k}}^{s}+n_{R,\vec{k}}^{s}n_{R,\vec{k}}^{s{\dagger}})+\tilde{m}_{N}(se^{i\theta_{\vec{k}}}n_{L,\vec{k}}^{s{\dagger}}n_{R,-\vec{k}}^{s{\dagger}}+se^{-i\theta_{\vec{k}}}n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s})]\,,$
(49)

which reproduces Eq.(42).

##  Appendix B Relevant time-evolution equations for two-point functions among
spinor fields and setting of those initial conditions

From Eq.(34) the time-evolution equations relevant to the BNA $n_{B}(t)$ in
Eq.(38) are derived as follows:

$\displaystyle{\cal D}_{t}\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (50)
$\displaystyle{\cal D}_{t}\langle n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
n_{R,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (51)
$\displaystyle{\cal D}_{t}\langle D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)$
$\displaystyle\hskip 14.22636pt-ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (52)
$\displaystyle{\cal D}_{t}\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{L,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle\right)$
$\displaystyle\hskip 14.22636pt-ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (53) $\displaystyle{\cal
D}_{t}\langle n_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+g\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (54)
$\displaystyle{\cal D}_{t}\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+g\langle
n_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt-ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{R,-\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,-\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (55)
$\displaystyle{\cal D}_{t}\langle
D_{R,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle+g\langle
n_{L,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (56)
$\displaystyle{\cal D}_{t}\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g\langle
n_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{R,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (57)
$\displaystyle{\cal D}_{t}\langle D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g\langle
n_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle\right)$
$\displaystyle\hskip 14.22636pt-
ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{L,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (58)
$\displaystyle{\cal D}_{t}\langle D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{R,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g\langle
n_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle\right)\,,$ (59)

$\displaystyle{\cal D}_{t}\langle D_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+g\langle
n_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (60)
$\displaystyle{\cal D}_{t}\langle D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{L,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)$
$\displaystyle\hskip 14.22636pt-ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g\langle
n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (61)
$\displaystyle{\cal D}_{t}\langle n_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
n_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{L,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
n_{L,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (62)
$\displaystyle{\cal D}_{t}\langle D_{R,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{R,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle+g\langle
D_{R,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle+g\langle
n_{L,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle\right)\,,$ (63)
$\displaystyle{\cal D}_{t}\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+g\langle
n_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt-ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
D_{L,-\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,-\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (64)
$\displaystyle{\cal D}_{t}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle+g\langle
n_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (65)
$\displaystyle{\cal D}_{t}\langle n_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
n_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle\right)$
$\displaystyle\hskip 14.22636pt-
ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{R,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g\langle
n_{R,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (66)
$\displaystyle{\cal D}_{t}\langle D_{L,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle$
$\displaystyle=2is\left|\vec{k}\right|\langle
D_{L,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle-
ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+g\langle
n_{R,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (67)
$\displaystyle{\cal D}_{t}\langle
D_{L,\vec{k}}^{s\dagger}D_{R,\vec{k}}^{s}\rangle$
$\displaystyle=ise^{-i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,-\vec{k}}^{s\dagger}\rangle+g\langle
D_{L,\vec{k}}^{s\dagger}n_{L,-\vec{k}}^{s\dagger}\rangle\right)$
$\displaystyle\hskip 14.22636pt-ise^{i\theta_{\vec{k}}}\left(m_{\chi}\langle
D_{R,-\vec{k}}^{s}D_{R,\vec{k}}^{s}\rangle+g\langle
n_{R,-\vec{k}}^{s}D_{R,\vec{k}}^{s}\rangle\right)\,,$ (68) $\displaystyle{\cal
D}_{t}\langle n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)$ $\displaystyle\hskip
14.22636pt+ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}\langle
n_{R,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g\langle
n_{R,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,.$ (69)

To solve the sets of equations derived as above, we also need to place the
initial conditions. We expect that the mass-eigenstate fields in the system
initially behave as free fields because the fields do not have any other
interactions except those due to $\langle\sigma(t)\rangle$. In that case, the
initial conditions on the mass-eigenstate fields are simply fixed through the
zero-particle states defined in the conventional static (adiabatic) vacuum.

We introduce a collective two-component spinor field as

$\left[\psi_{\vec{k}}^{s}\right]^{i}=\left(\begin{array}[]{c}n_{L,\vec{k}}^{s}(t)\\\
n_{R,\vec{k}}^{s}(t)\\\ D_{L,\vec{k}}^{s}(t)\\\
D_{R,\vec{k}}^{s}(t)\end{array}\right)^{i}\,.$ (70)

By defining

$\begin{aligned}
&{\left[A_{\vec{k}}\right]^{ij}\equiv\frac{1}{2V}\sum_{s}\left(\left\langle\left[\psi_{\vec{k}}^{s\dagger}\right]^{i}\left[\psi_{\vec{k}}^{s}\right]^{j}\right\rangle-\left\langle\left[\psi_{-\vec{k}}^{s}\right]^{j}\left[\psi_{-\vec{k}}^{s\dagger}\right]^{i}\right\rangle\right)}\\\
&{\left[\bar{A}_{\vec{k}}\right]^{ij}\equiv\frac{1}{2V}\sum_{s}s\left(\left\langle\left[\psi_{\vec{k}}^{s\dagger}\right]^{i}\left[\psi_{\vec{k}}^{s}\right]^{j}\right\rangle-\left\langle\left[\psi_{-\vec{k}}^{s}\right]^{j}\left[\psi_{-\vec{k}}^{s\dagger}\right]^{i}\right\rangle\right)}\\\
&{\left[B_{\vec{k}}\right]^{ij}\equiv\frac{1}{2V}\sum_{s}\left(\left\langle
se^{i\theta_{\vec{k}}}\left[\psi_{-\vec{k}}^{s}\right]^{i}\left[\psi_{\vec{k}}^{s}\right]^{j}\right\rangle+\left\langle
se^{i\theta_{\vec{k}}}\left[\psi_{-\vec{k}}^{s}\right]^{j}\left[\psi_{\vec{k}}^{s}\right]^{i}\right\rangle\right)}\\\
&{\left[\bar{B}_{\vec{k}}\right]^{ij}\equiv\frac{1}{2V}\sum_{s}s\left(\left\langle
se^{i\theta_{\vec{k}}}\left[\psi_{-\vec{k}}^{s}\right]^{i}\left[\psi_{\vec{k}}^{s}\right]^{j}\right\rangle+\left\langle
se^{i\theta_{\vec{k}}}\left[\psi_{-\vec{k}}^{s}\right]^{j}\left[\psi_{\vec{k}}^{s}\right]^{i}\right\rangle\right)}\end{aligned}\,,$
(71)

their initial conditions at $t=t_{0}$ are then set as Enomoto:2022gsx

$\begin{gathered}{\left[A_{\vec{k}}\left(t_{0}\right)\right]^{ij}=\left[\bar{B}_{\vec{k}}\left(t_{0}\right)\right]^{ij}=0}\,,\\\
{\left[\bar{A}_{\vec{k}}\left(t_{0}\right)\right]^{ij}=\sum_{\ell}\left[U^{*}\right]^{i\ell}\frac{|\vec{k}|}{\omega_{k}^{\ell}}\left[U^{T}\right]^{\ell
j}}\,,\\\
{\left[B_{\vec{k}}\left(t_{0}\right)\right]^{ij}=\sum_{\ell}[U]^{i\ell}\frac{\mathcal{M}_{d}^{\ell\ell}}{\omega_{k}^{\ell}}\left[U^{T}\right]^{\ell
j}}\,,\end{gathered}$ (72)

where $\mathcal{M}_{d}$ is a diagonalized mass matrix related to the mass
matrix $\mathcal{M}$ through a (orthogonal) transformation by $U$, which we
have performed the Takagi factorization method:

$\displaystyle\mathcal{M}$
$\displaystyle=\left(\begin{array}[]{cccc}0&\tilde{m}_{N}&0&g\\\
\tilde{m}_{N}&0&g&0\\\ 0&g&0&m_{\chi}\\\ g&0&m_{\chi}&0\end{array}\right)\,,$
(77) $\displaystyle\mathcal{M}_{d}$ $\displaystyle\equiv U^{T}\mathcal{M}U\,.$
(78)

In Eq.(72) $\omega_{k}^{\ell}$ is the energy with the diagonalized mass as

$\omega_{k}^{\ell}\equiv\sqrt{|\vec{k}|^{2}+\left(\mathcal{M}_{d}^{\ell\ell}\right)^{2}}.$
(79)

In terms of the two-point functions, the initial conditions (set at $t_{0}=0$)
read Enomoto:2022gsx

$\begin{array}[]{lll}\left\langle O_{1}\right\rangle_{t=0}&=\left\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(1+s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{11}\right)\,,\\\
\left\langle O_{2}\right\rangle_{t=0}&=\left\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{12}\,,\\\
\left\langle O_{3}\right\rangle_{t=0}&=\left\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{41}\,,\\\
\left\langle O_{4}\right\rangle_{t=0}&=\left\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{31}\right)\,,\\\
\left\langle O_{5}\right\rangle_{t=0}&=\left\langle
n_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(1+s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{22}\right)\,,\\\
\left\langle O_{6}\right\rangle_{t=0}&=\left\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{\boldsymbol{k}}\left(t_{0}\right)\right]^{42}\right)\,,\\\
\left\langle O_{7}\right\rangle_{t=0}&=\left\langle
D_{R,-\vec{k}}^{s+}D_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(1+s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{44}\right)\,,\\\
\left\langle O_{8}\right\rangle_{t=0}&=\left\langle
D_{R,\vec{k}}^{st}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{41}\right)\,,\\\
\left\langle O_{9}\right\rangle_{t=0}&=\left\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=-\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{32}\,,\\\
\left\langle O_{10}\right\rangle_{t=0}&=\left\langle
D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=-\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{42}\,,\\\
\left\langle O_{11}\right\rangle_{t=0}&=\left\langle
D_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{34}\,,\\\
\left\langle O_{12}\right\rangle_{t=0}&=\left\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{31}\,,\\\
\left\langle O_{13}\right\rangle_{t=0}&=\left\langle
n_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{11}\,,\\\
\left\langle O_{14}\right\rangle_{t=0}&=\left\langle
D_{R,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=-\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{44}\,,\\\
\left\langle O_{15}\right\rangle_{t=0}&=\left\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{32}\right)\,,\\\
\left\langle O_{16}\right\rangle_{t=0}&=\left\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(1+s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{33}\right)\,,\\\
\left\langle O_{17}\right\rangle_{t=0}&=\left\langle
n_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\right\rangle_{t=0}&=-\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{22}\,,\\\
\left\langle O_{18}\right\rangle_{t=0}&=\left\langle
D_{L,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{Vs}{2}\left[B_{k}\left(t_{0}\right)\right]^{33}\,,\\\
\left\langle O_{19}\right\rangle_{t=0}&=\left\langle
D_{L,\vec{k}}^{s\dagger}D_{R,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{34}\right)\,,\\\
\left\langle O_{20}\right\rangle_{t=0}&=\left\langle
n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle_{t=0}&=\frac{V}{2}\left(s\left[\bar{A}_{k}\left(t_{0}\right)\right]^{12}\right)^{*}\,.\end{array}$
(80)

## References

  * (1) P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016) doi:10.1051/0004-6361/201525830 [arXiv:1502.01589 [astro-ph.CO]].
  * (2) R. H. Cyburt, B. D. Fields, K. A. Olive and T. H. Yeh, Rev. Mod. Phys. 88, 015004 (2016) doi:10.1103/RevModPhys.88.015004 [arXiv:1505.01076 [astro-ph.CO]].
  * (3) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * (4) A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5, 32-35 (1967) doi:10.1070/PU1991v034n05ABEH002497
  * (5) A. D. Dolgov and D. P. Kirilova, Sov. J. Nucl. Phys. 51, 172-177 (1990) JINR-E2-89-321.
  * (6) J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 42, 2491-2504 (1990) doi:10.1103/PhysRevD.42.2491
  * (7) L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195-3198 (1994) doi:10.1103/PhysRevLett.73.3195 [arXiv:hep-th/9405187 [hep-th]].
  * (8) Y. Shtanov, J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 51, 5438-5455 (1995) doi:10.1103/PhysRevD.51.5438 [arXiv:hep-ph/9407247 [hep-ph]].
  * (9) L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 3258-3295 (1997) doi:10.1103/PhysRevD.56.3258 [arXiv:hep-ph/9704452 [hep-ph]].
  * (10) M. A. Amin, M. P. Hertzberg, D. I. Kaiser and J. Karouby, Int. J. Mod. Phys. D 24, 1530003 (2014) doi:10.1142/S0218271815300037 [arXiv:1410.3808 [hep-ph]].
  * (11) K. D. Lozanov, [arXiv:1907.04402 [astro-ph.CO]].
  * (12) A. G. Cohen and D. B. Kaplan, Phys. Lett. B 199, 251-258 (1987) doi:10.1016/0370-2693(87)91369-4
  * (13) A. Kusenko, K. Schmitz and T. T. Yanagida, Phys. Rev. Lett. 115, no.1, 011302 (2015) doi:10.1103/PhysRevLett.115.011302 [arXiv:1412.2043 [hep-ph]].
  * (14) M. Ibe and K. Kaneta, Phys. Rev. D 92, no.3, 035019 (2015) doi:10.1103/PhysRevD.92.035019 [arXiv:1504.04125 [hep-ph]].
  * (15) R. Daido, N. Kitajima and F. Takahashi, JCAP 07, 046 (2015) doi:10.1088/1475-7516/2015/07/046 [arXiv:1504.07917 [hep-ph]].
  * (16) F. Takahashi and M. Yamada, JCAP 10, 010 (2015) doi:10.1088/1475-7516/2015/10/010 [arXiv:1507.06387 [hep-ph]].
  * (17) P. Adshead and E. I. Sfakianakis, Phys. Rev. Lett. 116, no.9, 091301 (2016) doi:10.1103/PhysRevLett.116.091301 [arXiv:1508.00881 [hep-ph]].
  * (18) F. Takahashi and M. Yamada, Phys. Lett. B 756, 216-220 (2016) doi:10.1016/j.physletb.2016.03.020 [arXiv:1510.07822 [hep-ph]].
  * (19) A. Kusenko, L. Pearce and L. Yang, Phys. Rev. D 93, no.11, 115005 (2016) doi:10.1103/PhysRevD.93.115005 [arXiv:1604.02382 [hep-ph]].
  * (20) A. Maleknejad, JCAP 12, 027 (2016) doi:10.1088/1475-7516/2016/12/027 [arXiv:1604.06520 [hep-ph]].
  * (21) A. De Simone and T. Kobayashi, JCAP 08, 052 (2016) doi:10.1088/1475-7516/2016/08/052 [arXiv:1605.00670 [hep-ph]].
  * (22) A. De Simone and T. Kobayashi, JCAP 02, 036 (2017) doi:10.1088/1475-7516/2017/02/036 [arXiv:1610.05783 [hep-ph]].
  * (23) M. Son, F. Ye and T. You, Phys. Rev. D 99, no.9, 095016 (2019) doi:10.1103/PhysRevD.99.095016 [arXiv:1804.06599 [hep-ph]].
  * (24) A. Dasgupta, R. K. Jain and R. Rangarajan, Phys. Rev. D 98, no.8, 083527 (2018) doi:10.1103/PhysRevD.98.083527 [arXiv:1808.04027 [hep-ph]].
  * (25) K. J. Bae, J. Kost and C. S. Shin, Phys. Rev. D 99, no.4, 043502 (2019) doi:10.1103/PhysRevD.99.043502 [arXiv:1811.10655 [hep-ph]].
  * (26) V. Domcke, Y. Ema and K. Mukaida, JHEP 02, 055 (2020) doi:10.1007/JHEP02(2020)055 [arXiv:1910.01205 [hep-ph]].
  * (27) Y. P. Wu and K. Petraki, JCAP 01, 022 (2021) doi:10.1088/1475-7516/2021/01/022 [arXiv:2008.08549 [hep-ph]].
  * (28) M. Berbig, JHEP 01, 061 (2024) doi:10.1007/JHEP01(2024)061 [arXiv:2307.14121 [hep-ph]].
  * (29) A. Kusenko, L. Pearce and L. Yang, Phys. Rev. Lett. 114, no.6, 061302 (2015) doi:10.1103/PhysRevLett.114.061302 [arXiv:1410.0722 [hep-ph]].
  * (30) L. Yang, L. Pearce and A. Kusenko, Phys. Rev. D 92, no.4, 043506 (2015) doi:10.1103/PhysRevD.92.043506 [arXiv:1505.07912 [hep-ph]].
  * (31) L. Pearce, L. Yang, A. Kusenko and M. Peloso, Phys. Rev. D 92, no.2, 023509 (2015) doi:10.1103/PhysRevD.92.023509 [arXiv:1505.02461 [hep-ph]].
  * (32) H. Gertov, F. Sannino, L. Pearce and L. Yang, Phys. Rev. D 93, no.11, 115042 (2016) doi:10.1103/PhysRevD.93.115042 [arXiv:1601.07753 [hep-ph]].
  * (33) M. Kawasaki, L. Pearce, L. Yang and A. Kusenko, Phys. Rev. D 95, no.10, 103006 (2017) doi:10.1103/PhysRevD.95.103006 [arXiv:1701.02175 [hep-ph]].
  * (34) Y. P. Wu, L. Yang and A. Kusenko, JHEP 12, 088 (2019) doi:10.1007/JHEP12(2019)088 [arXiv:1905.10537 [hep-ph]].
  * (35) S. M. Lee, K. y. Oda and S. C. Park, JHEP 03, 083 (2021) doi:10.1007/JHEP03(2021)083 [arXiv:2010.07563 [hep-ph]].
  * (36) Y. Cado, B. von Harling, E. Massó and M. Quirós, JCAP 07, 049 (2021) doi:10.1088/1475-7516/2021/07/049 [arXiv:2102.13650 [hep-ph]].
  * (37) X. R. Wang, J. Y. Li, S. Enomoto, H. Ishida and S. Matsuzaki, Phys. Rev. D 108, no.2, 023512 (2023) doi:10.1103/PhysRevD.108.023512 [arXiv:2206.00519 [hep-ph]].
  * (38) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 06, 088 (2009) doi:10.1088/1126-6708/2009/06/088 [arXiv:0903.4155 [hep-lat]].
  * (39) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022) doi:10.1093/ptep/ptac097
  * (40) A. Deur, S. J. Brodsky and C. D. Roberts, Prog. Part. Nucl. Phys. 134, 104081 (2024) doi:10.1016/j.ppnp.2023.104081 [arXiv:2303.00723 [hep-ph]].
  * (41) R. Barate et al. [ALEPH], Z. Phys. C 76, 1-14 (1997) doi:10.1007/s002880050522
  * (42) C. X. Cui, J. Y. Li, H. Ishida, M. Kawaguchi, S. Matsuzaki and A. Tomiya, Universe 10 (2024) no.2, 68 doi:10.3390/universe10020068 [arXiv:2112.13533 [hep-ph]].
  * (43) R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338-341 (1984) doi:10.1103/PhysRevD.29.338
  * (44) G. Fejos and T. Hatsuda, [arXiv:2404.00554 [hep-ph]].
  * (45) L. Sagunski, P. Schicho and D. Schmitt, Phys. Rev. D 107, no.12, 123512 (2023) doi:10.1103/PhysRevD.107.123512 [arXiv:2303.02450 [hep-ph]].
  * (46) T. Bringmann, J. M. Cline and J. M. Cornell, Phys. Rev. D 99, no.3, 035024 (2019) doi:10.1103/PhysRevD.99.035024 [arXiv:1810.08215 [hep-ph]].
  * (47) B. Fornal and B. Grinstein, Phys. Rev. Lett. 120, no.19, 191801 (2018) [erratum: Phys. Rev. Lett. 124, no.21, 219901 (2020)] doi:10.1103/PhysRevLett.120.191801 [arXiv:1801.01124 [hep-ph]].
  * (48) J. M. Cline and J. M. Cornell, JHEP 07, 081 (2018) doi:10.1007/JHEP07(2018)081 [arXiv:1803.04961 [hep-ph]].
  * (49) B. Grinstein, C. Kouvaris and N. G. Nielsen, Phys. Rev. Lett. 123, no.9, 091601 (2019) doi:10.1103/PhysRevLett.123.091601 [arXiv:1811.06546 [hep-ph]].
  * (50) W. Y. Keung, D. Marfatia and P. Y. Tseng, JHEP 09, 053 (2019) doi:10.1007/JHEP09(2019)053 [arXiv:1905.03401 [hep-ph]].
  * (51) W. Y. Keung, D. Marfatia and P. Y. Tseng, JHEP 07, 181 (2020) doi:10.1007/JHEP07(2020)181 [arXiv:2001.09140 [hep-ph]].
  * (52) F. Elahi and M. Mohammadi Najafabadi, Phys. Rev. D 102, no.3, 035011 (2020) doi:10.1103/PhysRevD.102.035011 [arXiv:2005.00714 [hep-ph]].
  * (53) B. Fornal, B. Grinstein and Y. Zhao, Phys. Lett. B 811, 135869 (2020) doi:10.1016/j.physletb.2020.135869 [arXiv:2005.04240 [hep-ph]].
  * (54) B. Fornal and B. Grinstein, Mod. Phys. Lett. A 35, no.31, 2030019 (2020) doi:10.1142/S0217732320300190 [arXiv:2007.13931 [hep-ph]].
  * (55) S. Fajfer and D. Susič, Phys. Rev. D 103, no.5, 055012 (2021) doi:10.1103/PhysRevD.103.055012 [arXiv:2010.08367 [hep-ph]].
  * (56) G. Alonso-Álvarez, G. Elor, M. Escudero, B. Fornal, B. Grinstein and J. Martin Camalich, Phys. Rev. D 105, no.11, 115005 (2022) doi:10.1103/PhysRevD.105.115005 [arXiv:2111.12712 [hep-ph]].
  * (57) A. Strumia, JHEP 02, 067 (2022) doi:10.1007/JHEP02(2022)067 [arXiv:2112.09111 [hep-ph]].
  * (58) J. Berger and G. Elor, Phys. Rev. Lett. 132, no.8, 081002 (2024) doi:10.1103/PhysRevLett.132.081002 [arXiv:2301.04165 [hep-ph]].
  * (59) B. Fornal, Universe 9, no.10, 449 (2023) doi:10.3390/universe9100449 [arXiv:2306.11349 [hep-ph]].
  * (60) J. H. Liang, Y. Liao, X. D. Ma and H. L. Wang, JHEP 12, 172 (2023) doi:10.1007/JHEP12(2023)172 [arXiv:2309.12166 [hep-ph]].
  * (61) S. Gardner and M. Zakeri, Universe 10, no.2, 67 (2024) doi:10.3390/universe10020067 [arXiv:2311.13649 [hep-ph]].
  * (62) M. Bastero-Gil, T. Huertas-Roldan and D. Santos, [arXiv:2403.08666 [astro-ph.CO]].
  * (63) D. McKeen, M. Pospelov and N. Raj, Phys. Rev. D 103, no.11, 115002 (2021) doi:10.1103/PhysRevD.103.115002 [arXiv:2012.09865 [hep-ph]].
  * (64) I. Doršner, S. Fajfer, A. Greljo, J. F. Kamenik and N. Košnik, Phys. Rept. 641, 1-68 (2016) doi:10.1016/j.physrep.2016.06.001 [arXiv:1603.04993 [hep-ph]].
  * (65) X. R. Wang, J. Y. Li, S. Enomoto, H. Ishida and S. Matsuzaki, Phys. Rev. D 108, no.2, 023512 (2023) doi:10.1103/PhysRevD.108.023512 [arXiv:2206.00519 [hep-ph]].
  * (66) G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L8 (2023) doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  * (67) G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 952 (2023) no.2, L37 doi:10.3847/2041-8213/ace18b [arXiv:2306.16220 [astro-ph.HE]].
  * (68) J. Antoniadis et al. [EPTA], Astron. Astrophys. 678 (2023), A48 doi:10.1051/0004-6361/202346841 [arXiv:2306.16224 [astro-ph.HE]].
  * (69) J. Antoniadis et al. [EPTA and InPTA], Astron. Astrophys. 678 (2023), A49 doi:10.1051/0004-6361/202346842 [arXiv:2306.16225 [astro-ph.HE]].
  * (70) J. Antoniadis et al. [EPTA and InPTA:], Astron. Astrophys. 678 (2023), A50 doi:10.1051/0004-6361/202346844 [arXiv:2306.16214 [astro-ph.HE]].
  * (71) H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. Res. Astron. Astrophys. 23, no.7, 075024 (2023) doi:10.1088/1674-4527/acdfa5 [arXiv:2306.16216 [astro-ph.HE]].
  * (72) D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. Astrophys. J. Lett. 951, no.1, L6 (2023) doi:10.3847/2041-8213/acdd02 [arXiv:2306.16215 [astro-ph.HE]].
  * (73) D. J. Reardon, A. Zic, R. M. Shannon, V. Di Marco, G. B. Hobbs, A. Kapur, M. E. Lower, R. Mandow, H. Middleton and M. T. Miles, et al. Astrophys. J. Lett. 951 (2023) no.1, L7 doi:10.3847/2041-8213/acdd03 [arXiv:2306.16229 [astro-ph.HE]].
  * (74) P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * (75) C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No and K. Rummukainen, et al. JCAP 03, 024 (2020) doi:10.1088/1475-7516/2020/03/024 [arXiv:1910.13125 [astro-ph.CO]].
  * (76) V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435-2446 (2006) doi:10.1088/0264-9381/23/7/014 [arXiv:gr-qc/0512039 [gr-qc]].
  * (77) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, Class. Quant. Grav. 23, 4887-4894 (2006) [erratum: Class. Quant. Grav. 23, 7361 (2006)] doi:10.1088/0264-9381/23/15/008
  * (78) S. Kawamura, T. Nakamura, M. Ando, N. Seto, K. Tsubono, K. Numata, R. Takahashi, S. Nagano, T. Ishikawa and M. Musha, et al. Class. Quant. Grav. 23, S125-S132 (2006) doi:10.1088/0264-9381/23/8/S17
  * (79) K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011) [erratum: Phys. Rev. D 95, no.10, 109901 (2017)] doi:10.1103/PhysRevD.83.044011 [arXiv:1101.3940 [astro-ph.CO]].