††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
#  QCD preheating: New frontier of baryogenesis

Xin-Ru Wang  Center for Theoretical Physics and College of Physics, Jilin
University, Changchun, 130012, China  Jin-Yang Li  Center for Theoretical
Physics and College of Physics, Jilin University, Changchun, 130012, China
KEK Theory Center, IPNS, Tsukuba, Ibaraki 305-0801, Japan  Seishi Enomoto
School of Physics, Sun Yat-Sen University, Guangzhou 510275, China  Hiroyuki
Ishida  Center for Liberal Arts and Sciences, Toyama Prefectural University,
Toyama 939-0398, Japan  Shinya Matsuzaki  Center for Theoretical Physics and
College of Physics, Jilin University, Changchun, 130012, China

###### Abstract

We find that QCD can create the cosmological matter abundance via out-of-
equilibrium processes during the QCD phase transition, that is what we call
the QCD preheating, where the dynamic transition of the QCD vacuum
characterized by the quark condensate takes place instantaneously.  This
mechanism works when the Universe undergoes subsequent supercooled QCD
transition. We also find that the QCD preheating can work to create the baryon
asymmetry of the Universe if there is the new physics communicated with QCD.
These are new pictures of the thermal history around the QCD-phase transition
epoch, and thus the dynamic aspect of the QCD vacuum opens a new frontier to
explore low-scale matter generation such as baryogenesis. Pursuing the QCD
reheating era would also help deeply understanding the subatomic-scale physics
in the thermal history of the Universe.

##  I Introduction

The QCD phase transition is of importance to understand the origin of mass and
property for matter, i.e., nucleon, and has extensively been explored.  In
application to cosmology, however, QCD does not play major role to account for
cosmological matter abundances observed in the Universe today including the
baryon asymmetry of the Universe. This is simply because the QCD sector does
not have baryon number violation nor large enough CP violation, and will never
be out of thermal equilibrium at lower scales due to the strong coupling
nature. This is “folklore” that is widely accepted, and thus, people tend to
call for Beyond the Standard Model, which is usually setup irrespective to the
QCD phase transition #1#1#1 QCD baryogenesis has been argued in Ipek:2018lhm ;
Servant:2014bla ; Croon:2019ugf  with higher scale confinement, not at the
ordinary QCD scale..

In this paper, we claim that dynamic aspects of the QCD vacuum during the QCD
phase transition can be crucial to create cosmological matter abundances
observed in the Universe today: the time-varying light-quark condensate
$\langle\bar{q}q\rangle$ produces the total matter-antimatter abundance,
through undergoing a supercooling and reheating before the static hadron phase
is formed.

The dynamic aspect of vacuum in quantum field theory has extensively been
explored so far in light of the reheating process after the inflation in the
early 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 ; Enomoto:2020aapps ). The
oscillating background field of the inflaton induces the nonadiabatic state
for other coupled species. It drives the nonperturbative particle production,
and the amount is exponentially amplified by the parametric resonance similar
to a swing with a pumping oscillator.

The “ballpark” of preheating has also been extended to the baryogenesis via
the nonperturbative production at high scales Funakubo:2000us ;
Rangarajan:2001yu ; Pearce:2015nga ; Enomoto:2017rvc ; Wu:2019ohx ;
Enomoto:2020lpf ; Lee:2020yaj . The mechanism of preheating is applied to
interesting variant scenarios: the production of massive particles heavier
than the inflaton Kolb:1996jt ; Kolb:1998he ; Giudice:1999fb ; Peloso:2000hy ;
the preheating due to the alternative field instead of the inflaton Garcia-
Bellido:2008ycs ; Adshead:2015kza ; Ema:2016dny ; Adshead:2018oaa , etc.

Preheating at the subatomic scale has never been discussed in the context of
baryogenesis. A scalar condensate is present even in QCD, that is the light
quark condensate $\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
symmetry and its breaking by the current quark masses. Therefore, the dynamic
motion of $\langle\bar{q}q\rangle$ should have the potential to explosively
produce the number densities for nucleon and anti-nucleon by nonadiabatic
processes, similarly to the preheating induced by the nonadiabatic-varying
vacuum. In this paper, we call the observation above the QCD preheating and
find it indeed work.

Remarkably enough, the QCD preheating reveals a new picture of the thermal
history related to nucleon: cosmological matters are produced out of
equilibrium and overwhelm the thermal abundance, which survives long enough
during the reheating epoch. This new picture of the thermal history is
applicable to the baryogenesis at the subatomic scale once the QCD preheating
is coupled to Beyond the Standard Model, that is the QCD-scale baryogenesis.

##  II Overview of the QCD preheating scenario

Figure 1: A schematic view of the QCD preheating in the thermal history of the
Universe. The top panel depicts the evolution of the temperature of the
Universe at around the QCD phase transition, where the QCD preheating via the
dynamic $\langle\sigma\rangle\sim\langle\bar{q}q\rangle$ happens at, say,
$T=T_{n}\sim 0.7T_{\rm pc}\simeq 109$ MeV by starting the roll-down in the
potential ($V$) passing the supercooling, and reheats the Universe up to
$T=T_{\rm reh}\simeq 119$ MeV (see the text).  The bottom panel traces the
time evolution of $T$ in terms of the potential shapes of the dynamic
$\langle\bar{q}q\rangle$ observed in the reference frame comoving with the
Hubble parameter. Since the $\langle\sigma\rangle$ motion is much faster than
the Hubble expansion rate, a sudden discontinuous jump-up of
$\langle\sigma\rangle$, at $T=T_{n}$ would be detected by eyes of observers
getting on the Hubble reference frame, which can be interpreted as a first
order transition in the thermal history of the Universe.  Passing the
relaxation as in Eq.(10), the system comes back to the normal thermal
equilibrium, which makes the potential of $\langle\sigma\rangle$ lifted up in
the standard manner of thermal QCD.

To begin with, we outline the static aspect of the QCD chiral phase
transition, and  turn to a possible dynamic feature of the QCD vacuum.

The lattice QCD has confirmed that thermal QCD with 2 + 1 flavors at the
physical point undergoes the crossover for the chiral phase transition, what
is called the chiral crossover, at the pseudocritical temperature $T_{\rm
pc}\sim 155$ MeV Aoki:2009sc ; Borsanyi:2011bn ; Ding:2015ona ;
Bazavov:2018mes ; Ding:2020rtq . At almost the same temperature, the
deconfinement-confinement transition (crossover) is expected to happen as well
Bazavov:2016uvm ; Ding:2017giu . Above $T_{\rm pc}$, the light quark
condensate $\langle\bar{q}q\rangle$ takes nearly vanishing values when
properly renormalized to be divergent free. Cooling down below $T_{\rm pc}$,
$\langle\bar{q}q\rangle$ starts to get sizable and finally reaches the value
measured at the vacuum when $T$ goes well below $T_{\rm pc}$, $T/T_{\rm
pc}\lesssim 0.7$ Aoki:2009sc .

This is the static picture of the QCD phase transition at our best knowledge.
When Hubble evolutionary Universe and some possible Beyond the Standard Model
contributions are taken into account, the dynamic picture of the QCD phase
transition might drastically be altered. In particular, subsequent
supercooling including the QCD phase transition can be realized when the
Standard Model is extended to be scale invariant along with a dark sector, as
discussed in Pisarski:1983ms ; Iso:2017uuu ; Hambye:2018qjv ; Sagunski:2023ynd
#2#2#2 In the literature Kim:2018byy ; Kim:2018knr ; Kim:2019lks , it was
argued that even without Beyond the Standard Model, a barrier-less boundary
for a single bubble between the chiral symmetric or deconfining and broken or
confining phases is interpreted as the crossover, which starts to happen at
$T=T_{\rm pc}$, that is consistent with what the lattice simulations have
observed. The confining hadron-phase bubbles expand until all bubbles occupy
the Universe, at which time, the so-called the nucleation and percolation,
$\langle\bar{q}q\rangle$ reaches the true vacuum (at $T\lesssim 0.7T_{\rm
pc}$). This implies that the evolution of the hadron phase acts as if the
Universe undergoes a supercooling in the regime $0.7T_{\rm pc}\lesssim
T\lesssim T_{\rm pc}$. The temperature where the hadron phase bubbles fully
cover the Universe, $T_{f}$, has been estimated to be indeed lower than
$T_{\rm pc}$: $T_{f}/T_{\rm pc}\sim 0.76$. However, in the literature
Kim:2018byy ; Kim:2018knr ; Kim:2019lks  the used pressure expressions are not
widely accepted, hence do not match with the lattice QCD thermodynamics. .

In the literature Sagunski:2023ynd  the QCD supercooling has been realized by
coupling to a dark sector which undergoes an ultrasupercooling, analogously to
the scenario of the QCD-induced supercooling electroweak phase transition
Iso:2017uuu . This subsequent supercooling is triggered by the highly
dominated potential energy density of the dark sector in the Hubble parameter,
which can be operative until $T\sim T_{\rm pc}$. All the massless six quark
condensates are trapped by the thermal potential at the origin (the chiral
symmetric phase). After the supercooling ends in the dark and the Higgs
sectors (via tunneling into the true vacua), quark mass terms are generated as
usual, and only lightest flavors (say, up and down quarks) survive in the QCD
thermal plasma. Thus the QCD sector also exits from the supercooling.

Then the light-quark condensate $\sigma\sim\langle\bar{q}q\rangle$ starts to
roll down to the true vacuum, along the generated linear (negative) slope
$\sim-m_{q}\sigma$ at the origin. This rolling motion is much faster than the
Hubble rate, so can be nonadiabatic and out-of-equilibrium to provide the
chance of particle production in the QCD phase transition epoch. This is the
setup that we are concerned about in the present paper.

Realization of such subsequent supercooling involves a nonperturbative
analysis on QCD with massless six flavors. In the literature Sagunski:2023ynd
, the authors have worked on chiral effective models in the mean field
approximation (only including one-loop corrections from quarks), which
corresponds to the limit of large number of colors (so-called the large
$N_{c}$ limit) and assumed that the usual lightest three flavors predominantly
form the quark condensates faster than other quarks. Then the nucleation
temperature $T_{n}$ for the completion of the QCD phase transition over the
Hubble evolution has been estimated as $T_{n}/T_{c}\sim 0.98$ for a
phenomenologically viable dark sector setup, where $T_{c}$ denotes the
critical temperature for the first order chiral phase transition in QCD with
massless three flavors. This value could be corrected by amount of ${\cal
O}(1/N_{c})\sim 30\%$ when the subleading contributions in the large $N_{c}$
expansion could be incorporated, which is to be $T_{n}/T_{c}\sim 0.7$.
Furthermore, a larger number of flavors could make the strength of the first
order phase transition stronger, which would imply the smaller $T_{n}/T_{c}$.
Thus we may suspect that the QCD supercooling could last long enough and
$T_{n}$ is low enough that the thermal correction to $\langle\bar{q}q\rangle$
can safely be neglected when it starts to roll down.

In this paper, we discuss the dynamic $\langle\bar{q}q\rangle$ by assuming a
possible supercooling-like picture during the QCD phase transition epoch as
above driven by coupling to the Beyond the Standard Model. The dynamic motion
of $\langle\bar{q}q\rangle$ instantaneously starts at say, $T\sim 0.7T_{\rm
pc}$ passing the supercooling. Since the scalar field $\langle\bar{q}q\rangle$
couples to the meson and baryon states, the nonadiabatic particle production
will then take place and reheat the Universe from $0.7T_{\rm pc}$ up to
$T_{\rm reh}(<T_{\rm pc})$. This is how the QCD preheating works and can
successfully create various matters and antimatters through the
nonperturbative and nonadiabatic processes, as will be explicitly seen below.

A schematic view of the QCD preheating is illustrated in Fig. 1.

##  III Dynamic QCD vacuum via a chiral effective theory description

To describe the dynamic $\langle\bar{q}q\rangle$ after the QCD supercooling
ends as elaborated in Sec.II, we employ a chiral effective theory, called the
linear sigma model, with the lightest two flavors, in which color-singlet
composite operators $\bar{q}q$ and $qqq$ are monitored by hadronic-
interpolating fields. The model is built based on the chiral $SU(2)_{L}\times
SU(2)_{R}$ symmetry and its breaking structure with the current quark mass for
the up and down quarks included. The building blocks are: i) the 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$); ii)
the 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 linear-
sigma model Lagrangian is thus given as

$\displaystyle\mathcal{L}=\operatorname{tr}\left[\partial_{\mu}M^{\dagger}\partial^{\mu}M\right]-V$
$\displaystyle\hskip 20.0pt+\bar{N}i\partial
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)

$f_{\pi}$ is the pion decay constant $\simeq 92.4$ MeV, $m_{\pi}$ the pion
mass $\simeq 140$ MeV, and $m_{N}$ the nucleon mass $\simeq 940$ MeV. The
vacuum expectation value of $\sigma$, $\langle\sigma\rangle$, and its
dynamical evolution are identical to those for $\langle\bar{q}q\rangle$. We
fix the potential parameters $m^{2}$ and $\lambda$ to the values determined at
the true vacuum $\langle\sigma\rangle=f_{\pi}$ satisfying the stationary
condition. Then we have $\lambda f_{\pi}^{2}=(M_{\sigma}^{2}-m_{\pi}^{2})/2$
and $m^{2}=3m_{\pi}^{2}/2-M_{\sigma}^{2}/2$, with the $\sigma$ mass squared
defined as
$M_{\sigma}^{2}=\partial^{2}V/\partial\sigma^{2}|_{\sigma=f_{\pi}}$, which we
take $\simeq(500\>{\rm MeV})^{2}$.

The equation of motion for $\langle\sigma\rangle\sim\langle\bar{q}q\rangle$
with the space homogeneity leads

$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\,.$
(2)

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. $\gamma$ plays the role of the full width of the $\sigma$
meson identified as $f_{0}(500)$ in the Particle Data Group. 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:2020ssz .

The dynamic oscillation of $\langle\bar{q}q\rangle$ will globally be processed
along the potential, which is well approximated to be $V$ in Eq.(1) after the
supercooling ends at, say, $T_{n}=0.7T_{\rm pc}$. Since the oscillator
$\langle\bar{q}q\rangle$ couples to the meson and baryon states, the
nonadiabatic particle production will then take place and reheat the universe
from $T=T_{n}$ up to the reheating temperature $T_{\rm reh}(<T_{\rm pc})$. The
typical time scale (1/500, MeV-1) of the dynamic $\langle\bar{q}q\rangle$
oscillation is much faster than the background Hubble evolution ($t\sim
1/H\sim 1/(T^{2}/M_{p})$ with $T\sim 0.1\,{\rm GeV}$ and $M_{p}\sim
10^{19}\,{\rm GeV}$ being the Planck scale). Hence the Hubble Universe will
observe as if $\langle\bar{q}q\rangle$ instantaneously jumps to the true
vacuum (hadron phase) from the origin (symmetric phase).

The backreaction terms correspond to the effective mass of $\sigma$ arising as
the plasma effect through the nonadiabatic production of the pion and the
nucleons due to the dynamic $\sigma$. We have checked that their effect on the
motion of $\langle\sigma\rangle$ is sufficiently small. The nonadiabatic
particle production occurs around $\langle\sigma\rangle\sim 0$ at once, but
the parametric resonance does not since the later motion does not come back
there due to the friction $\gamma$. Thus the dynamic $\langle\sigma\rangle$ is
well approximately described only by Eq.(2) without interaction terms with
pions and nucleons. The time evolution of $\langle\sigma\rangle$ is plotted in
Fig. 2, as a function of $(M_{\sigma}t)$.

Figure 2: The time evolution of $\langle\sigma\rangle$ (blue) and the total
nucleon number density $(n_{N}+\bar{n}_{N})$ (red), with the initial condition
$\langle\sigma(0)\rangle=\langle\dot{\sigma}(0)\rangle=0$, inferred from the
QCD preheating, and $n_{N}(0)=\bar{n}_{N}(0)=0$. $\langle\sigma\rangle$
arrives at the true vacuum, $\langle\sigma\rangle=f_{\pi}=92.4$ MeV, at
$M_{\sigma}t\sim 8$. The nonadiabatic production of $(n_{N}+\bar{n}_{N})$ ends
when $M_{\sigma}t\sim 5$.

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

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

where we have assumed the relativistic degrees of freedom $g_{*}\simeq 14$ at
$T=0.7T_{\rm pc}\simeq 109\>{\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(119\>{\rm MeV})^{3}\simeq
1.24\times 10^{7}\>{\rm MeV}^{3}.$ (4)

##  IV Nonadiabatic production of nucleons due to the dynamic
$\langle\sigma\rangle$

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.$
(5)

Thus 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$ depicted in Fig.
2. Such a time-varying mass 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},$ (6)

where we read
$\langle\dot{\sigma}\rangle\simeq\sqrt{V(\sigma=0)-V(\sigma=f_{\pi})}\sim(135\>{\rm
MeV})^{2}$. Hence the nonperturbative nucleon production would be completed
within the range in Eq. (6). In terms of Fig. 2, it corresponds to

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

The actual production time can be earlier than the number in Eq.(7) because
the estimated velocity $\langle\dot{\sigma}\rangle$ would be smaller due to
the friction $\gamma$.

The total nucleon number density can be evaluated as #3#3#3 $\vec{A}$ denotes
the spatial component of the vector variable $A$.

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

for $N=p,n$. Here $\tilde{\rho}_{N}(\vec{k},t)$ is the kinetic energy density
of the nucleon in momentum space, which is derived from the Hamiltonian as

$\displaystyle\tilde{\rho}_{N}(\vec{k},t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}\left[\bar{N}_{L}(\vec{k},L)\left(\vec{\gamma}\cdot\vec{k}+\tilde{m}_{N}(t)\right)N_{L}(\vec{k},t)\right.$
(9)
$\displaystyle\quad\left.+\bar{N}_{R}(\vec{k},t)\left(\vec{\gamma}\cdot\vec{k}+\tilde{m}_{N}(t)\right)N_{R}(\vec{k},t)\right]$
$\displaystyle+2\omega(\vec{k},t)\,,$

with $\omega(k,t)=\sqrt{|\vec{k}|^{2}+\tilde{m}_{N}^{2}(t)}$ being the one-
particle energy of the nucleon, $V=\int d^{3}x$ the space volume of the
system, and $N_{L}(\vec{k},t)$ and $N_{R}(\vec{k},t)$ the Fourier transformed
Dirac fields. The last term in Eq.(9) 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.(8) is thus determined by
solving coupled equations of motion for $N_{L}(\vec{k},t)$,
$N_{R}(\vec{k},t)$, and $\langle\sigma(t)\rangle$, which is plotted in Fig. 2.
For more details, see Appendix A.

We find that the total number density $(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$.  Note that this number density is much larger
than the thermal equilibrium density at the reheating temperature $T_{\rm
reh}\simeq 119$ MeV, $[(n_{N}+\bar{n}_{N})]_{\rm EQ}\sim 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
$M_{\sigma}\Delta t_{\rm relax}=M_{\sigma}(n_{N}\langle\sigma
v\rangle)^{-1}\sim 700$, or equivalently,

$\displaystyle\Delta t_{\rm relax}\sim 5\times 10^{-7}\ {\rm fs}\,,$ (10)

where 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 stock 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.

This fact tempts us to consider the application to the baryogenesis, in which
the Sakharov criteria are required: the baryon number violating interaction,
C- and CP-violating interactions, and out-of-equilibrium condition
Sakharov:1967dj .

Once the baryonic asymmetry
$\epsilon\equiv(n_{N}-\bar{n}_{N})/(n_{N}+\bar{n}_{N})$ is provided by the
other mechanism until the relaxation time, the QCD preheating can successfully
generate the desired amount of the net baryon number of the Universe,

$\displaystyle Y_{B}\equiv\frac{(n_{N}+\bar{n}_{N})}{s}\cdot\epsilon\sim
10^{-10}\times\left(\frac{\epsilon}{10^{-8}}\right)\,,$ (11)

where the entropy density in Eq.(4) has been taken into account.

Here a salient feature is seen: the observed baryon number can be realized by
a relatively smaller asymmetry than that accumulated in the thermal
equilibrium due to the overproduced nucleons, which will not be washed out as
long as the asymmetry generation completes before the system comes back to the
thermal equilibrium.

It is true that QCD in the Standard Model cannot solely generate the asymmetry
$\epsilon$ in Eq.(11) because of the absence of sufficient CP- and baryon
number-violating interactions, but Eq.(11) opens a new roadmap as new
baryogenesis involving a number of the new sector candidates coupled to the
Standard Model: once an external sector transferring the CP and baryon number
violations to QCD is hypothesized, the QCD preheating triggered by the dynamic
$\langle\bar{q}q\rangle$ can create the net baryon asymmetry through Eq.(11),
if and only if it ends until the relaxation time $M_{\sigma}\Delta t_{\rm
relax}\sim 700$.

##  V Production of baryon asymmetry

As a benchmark, we introduce a class of models in which the baryogenesis makes
the most of the QCD preheating to supply $\epsilon$ in Eq.(11).

Consider a dark sector with dark-Majorana neutron fields $n_{D_{L,R}}$ being
allowed to couple to the QCD neutron fields $n_{L,R}$ in a minimal way #4#4#4
The mixing angle between the dark sector and QCD via the $g_{L,R}$ couplings
are of ${\cal O}(10^{-4})$ for benchmark models in Fig. 3, so that the dark
sector contributions to the successful QCD hadron physics and chiral phase
transition are safely negligible. :

$\displaystyle{\cal L}_{n\mathchar 32\relax n_{D}}$
$\displaystyle=-m_{D}\bar{n}_{D_{R}}n_{D_{L}}-\frac{1}{2}M_{L}\overline{n_{D_{L}}^{c}}n_{D_{L}}-\frac{1}{2}M_{R}\overline{n_{D_{R}}^{c}}n_{D_{R}}$
$\displaystyle\quad-g_{L}\bar{n}_{R}n_{D_{L}}-g_{R}\bar{n}_{D_{R}}n_{L}+\text{
h.c.}\,,$ (12)

where the superscript $c$ stands for the charge conjugation. In general, only
two CP-violating phases can be physical to be introduced among the mass
couplings ($m_{D},g_{L},g_{R},M_{L},M_{R}$).

In the dark sector scenario described in Eq.(12), the dynamic
$\langle\sigma(t)\rangle$ also plays an important role in generating the CP-
violating source. Indeed, even the physical CP phases on the mass parameters
can be erased by the appropriate diagonalization of the mass matrix since the
Lagrangian consists of only the two-point interactions. However, the neutron
mass currently depends on time via $\langle\sigma(t)\rangle$ as in Eq.(5), and
thus the CP phases can re-appear from the kinetic terms as the so-called Berry
connection. A similar discussion can be seen in Enomoto:2018yeu . The combined
effect of the CP phase and the baryon number violating couplings thus enables
the simultaneous asymmetric production of the neutron and antineutron by the
dynamic $\langle\sigma\rangle$.

For simplicity, we turn off one of two phases and embed it in $M_{L}$. With
this phase, the nonzero Majorana mass couplings thus transfer the CP violation
and the baryon number violation into the QCD sector. The dark baryons need to
be as heavy as the neutron, in such a way to make the dark-sector
communication with the neutron operative in the particle production process.

A similar idea to transfer the baryon number violation from a dark sector into
the QCD sector via introducing a neutron-dark neutron coupling has been
discussed in the literature Bringmann:2018sbs , where the production mechanism
of the number densities differs from the nonperturbative one, however.

We shall evaluate the net baryon number instead of the asymmetry $\epsilon$,
which can be defined by the (approximately conserved) $U(1)$ Noether charge as
#5#5#5 In the present dynamic system, it is actually less costly to work on
the evaluation of the dynamical net nucleon number instead of a direct
estimate of $\epsilon$.

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

where the last term corresponds to the subtraction of the divergent part
induced by the zero-point energy ($4\times\frac{1}{2}$). The time evolution of
the vacuum expectation values in Eq.(13) is evaluated by solving chained
equations of motion for $n_{L,R}$ coupled to $n_{D_{L,R}}$ through Eq.(12),
together with the dynamic $\langle\sigma(t)\rangle$ obeying Eq.(2).

The definition of the net baryon number in Eq.(13) should potentially include
the proton contribution as well, which, however, we can safely ignore in our
analysis. This is because the proton does not have the CP- and baryon number-
violating interactions, and thus the net baryon number by the proton cannot be
generated.

The sigma motion makes the net baryon number produced through the baryon-
number violating couplings $g_{L,R}$. We have observed that the net number
starts to oscillate even after the baryon number production. This is because
of the presence of the $n$-$n_{D}$ and $n$-$\bar{n}$ oscillations that last
eternally. To obtain the static net baryon number, the couplings $g_{L,R}$
connecting the neutrons and the dark sector need to somehow get damping in
time or vanishing at the later era, such as $g_{L,R}(t)=\varphi(t)g_{L0,R0}$
with $\varphi(t)=e^{-\Gamma_{\varphi}t}\cos m_{\varphi}t$. In that case, since
the magnitude of the $g_{L,R}$ couplings asymptotically and promptly drops to
zero, any phenomenological and astrophysical bounds can safely be satisfied
when the observation is made.

The damping-oscillation factor $\varphi(t)$ could arise when one considers an
underlying picture of the $g_{L,R}$ mass-mixing coupling to be given by a
scalar or dilaton, which communicates between the dark sector and the normal
QCD sector. In that case, $\varphi(t)$ would be regarded as a background part
of the scalar field, and the form of the damping oscillation in time would be
provided by the scalar decay width $\Gamma_{\varphi}$, following the time
evolution equation:
$\langle\ddot{\varphi}\rangle+\Gamma_{\varphi}\langle\dot{\varphi}\rangle+m_{\varphi}^{2}\langle\varphi\rangle=0$.
This ultraviolet completion could be rich in phenomenology and cosmology,
which is, however, beyond the current scope, to be pursued elsewhere.

Figure 3 shows the time evolution of the created net-baryon number yield
$Y_{B}=n_{B}/s$, normalized by the entropy density $s=10^{7}\,{\rm MeV}^{3}$
[see Eq.(4)]. The plot has been made for $M_{R}=1000\,{\rm MeV}$,
$M_{L}=1000\cdot e^{i\frac{\pi}{3}}\,{\rm MeV}$, $m_{D}=1050$ MeV,
$m_{N}=1000$ MeV, $m_{\varphi}=10$ MeV, and $\Gamma_{\varphi}=5$ MeV, along
with various sets of the initial $g_{L,R}$ coupling strengths,
$(g_{L0},g_{R0})/{\rm MeV}=(1,2),(5,10)$, and $(10,20)$. The detailed tools
and methods for the present numerical computation are supplemented in Appendix
A.

We have found that the net baryon number scales like
$n_{B}\propto|g_{L}g_{R}|^{2}$. The net baryon number is essentially produced
by the nonadiabatic $\sigma$ motion and the associated nonperturbative
nucleon-antinucleon pair production for the time scale in the first phase
$0\leq M_{\sigma}t\lesssim 5$, which happens before $\langle\sigma\rangle$
reaches the true vacuum at $M_{\sigma}t\sim 10$ [See Fig. 2]. This production-
time range is consistent with a rough estimation done in Eq.(7). The time
evolution of $Y_{B}$ turns to the second phase ($5\lesssim M_{\sigma}t\lesssim
100$), where the intrinsic neutron-antineutron oscillation governs the
dynamics. In the end, coming to the third phase ($M_{\sigma}t\gtrsim 100$),
$Y_{B}$ asymptotically gets saturated to the constant due to the damping
oscillation of the couplings $g_{L,R}$, which start to operate at around
$M_{\sigma}t\sim M_{\sigma}/\Gamma_{\varphi}=100$. The general feature of
$Y_{B}$ as viewed in Fig. 3 is qualitatively independent of $m_{\varphi}$ when
it is lower compared with the particle production time scale $M_{\sigma}/5$
because $g_{L}$ and $g_{R}$ can maintain almost constants during the net
baryon number production. On the other hand, the higher mass, $m_{\varphi}\gg
M_{\sigma}/5$, might spoil the net number production because $g_{L}$ and
$g_{R}$ become ineffective due to their violent vibrations around zero.

The asymmetry production is thus completed faster than the aforementioned
relaxation time ($M_{\sigma}\Delta t_{\rm relax}\sim 700$). Thereby the
baryogenesis is successfully accomplished to yield the desired amount of the
baryon asymmetry, $Y_{B}=n_{B}/s\sim 10^{-10}-10^{-9}$. We would remark that
the typical baryogenesis scenario requires a high-scale physics, such as much
higher than the electroweak scale, whereas the present one is realized at the
lower scale due to the QCD preheating.

Figure 3: The time evolution of the baryon asymmetry yield $Y_{B}=n_{B}/s$,
normalized with the entropy density $s=10^{7}\,{\rm MeV}^{3}$. The model
parameters have been fixed as described in the text.

##  VI Conclusions

The QCD preheating provides a nonperturbative cosmological particle production
out of equilibrium and reheats the Universe (Fig. 1). Consequently, the baryon
asymmetry of the Universe can successfully be yielded through Eq.(11), once a
source of the asymmetry is provided before the system comes back to the
thermal equilibrium, which corresponds to the time scale $\Delta t_{\rm
relax}\sim 10^{-7}\ {\rm fs}$ after the QCD preheating is completed.

The QCD preheating makes it possible to create not only the ordinary matters,
but also even the net baryon number or non-standard particles if there is the
new physics communicated with QCD.

The presently proposed QCD-scale baryogenesis can satisfy one of the most
difficult Sakharov criteria, namely out-of-equilibrium condition, without
introducing new physics effects. This is in sharp contrast to the conventional
way invoking the first-order electroweak phase transition or heavy particle
decays. Even the baryon number violation can be achieved around the QCD scale,
not higher scales. Thus the QCD preheating does substitute the conventional
methodologies to create the cosmological matter abundance observed in the
present-day Universe, including the baryon asymmetry of the Universe, hence is
a new paradigm for baryogenesis and paves the way for a new frontier involving
variant Beyond the Standard Model candidates.

Another QCD baryogenesis has been addressed based on higher scale features of
QCD, not at the subatomic scale Ipek:2018lhm ; Servant:2014bla ; Croon:2019ugf
, where the production mechanism of the number densities is neither due to the
nonperturbative time-dependent varying vacuum.

We introduced a class of the benchmark scenario to incorporate sufficient C-
and CP-violating and baryon number violating interactions just as a reference,
in which the generation of the baryonic asymmetry makes the most of the QCD
preheating, and showed that the baryogenesis works (Fig. 3).  Full embedding
of the dark and the Standard Model sectors into the scalegenesis would be
worth pursuing, which is to be explored elsewhere.

The QCD preheating leads to a sort of extra reheating in the thermal history,
as illustrated in Fig. 1, hence has an impact also on exploring thermal or
cosmological scenarios addressing the epoch around the QCD phase transition.
It would be, for instance, worth revisiting the estimate of the cosmological
abundance of the QCD axion, along with the QCD preheating.

The QCD preheating follows presence of the QCD supercooling as described in
Sec.II, which would be subsequent supercooling involving the Standard Model
Higgs and dark sectors with a scale invariant setup. Precise realization of
the QCD supercooling thus requires nonperturbative analysis on the vacuum
structure in massless-six flavor QCD and its tunneling rate. This issue
deserves to another publication.

## Acknowledgements

We thank Mamiya Kawaguchi for useful comment on the QCD phase transition. This
work was supported in part by the National Science Foundation of China (NSFC)
under Grant No.11975108, 12047569, 12147217 and the Seeds Funding of Jilin
University (S.M.).

##  Appendix A Preliminaries for numerical evaluation of time evolution of the
baryon number

In this Appendix we give a list of tools for the analysis on the QCD
preheating on the basis of the linear sigma model description and a dark
sector introduced in the main text. Note that we formulate the fermion sector
by the two-component spinors. The relations to the Dirac fermions are
followings;

$n=\left(\begin{array}[]{c}n_{L}\\\ n_{R}^{c\dagger}\end{array}\right),\qquad
n_{D}=\left(\begin{array}[]{c}n_{D_{L}}\\\
n_{D_{R}}^{c\dagger}\end{array}\right),$ (14)

where the superscript “$c$” denotes the charge conjugate.

###  A.1 Equations of motion

The relevant equations of motion are:

$\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}n_{L}-\tilde{m}_{N}^{*}n_{R}^{c\dagger}-g_{R}n_{D_{R}}^{c\dagger}\,,$
(15) $\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}n_{L}^{\dagger}-\tilde{m}_{N}n_{R}^{c}-g_{R}^{*}n_{D_{R}}^{c}\,,$
$\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}n_{R}^{c}-\tilde{m}_{N}^{*}n_{L}^{\dagger}-g_{L}^{*}n_{D_{L}}^{\dagger}\,,$
$\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}n_{L}^{\dagger}-\tilde{m}_{N}n_{R}^{c}-g_{L}n_{D_{L}}\,,$
$\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}n_{D_{L}}-m_{D}^{*}n_{D_{R}}^{c\dagger}-M_{L}^{*}n_{D_{L}}^{\dagger}-g_{L}^{*}n_{R}^{c\dagger}\,,$
$\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}n_{D_{L}}^{\dagger}-m_{D}n_{D_{R}}^{c}-M_{L}n_{D_{L}}-g_{L}n_{R}^{c}\,,$
$\displaystyle 0$
$\displaystyle=\bar{\sigma}^{\mu}i\partial_{\mu}n_{D_{R}}^{c}-m_{D}^{*}n_{D_{L}}^{\dagger}-M_{R}^{*}n_{D_{R}}^{c\dagger}-g_{R}n_{L}^{\dagger}\,,$
$\displaystyle 0$
$\displaystyle=\sigma^{\mu}i\partial_{\mu}n_{D_{R}}^{c\dagger}-m_{D}n_{D_{L}}-M_{R}n_{D_{R}}^{c}-g_{R}^{*}n_{L}\,,$

where $\tilde{m}_{N}(t)=m_{N}\frac{\langle\sigma(t)\rangle}{f_{\pi}}$.
Equivalently, the above equations can collectively be represented by

$\bar{\sigma}^{\mu}i\partial_{\mu}\left(\begin{array}[]{c}n_{L}\\\
n_{R}^{c}\\\ n_{D_{L}}\\\
n_{D_{R}}^{c}\end{array}\right)=\mathcal{M}^{*}\left(\begin{array}[]{c}n_{L}^{\dagger}\\\
n_{R}^{c\dagger}\\\ n_{D_{L}}^{\dagger}\\\
n_{D_{R}}^{c\dagger}\end{array}\right),\qquad\sigma^{\mu}i\partial_{\mu}\left(\begin{array}[]{c}n_{L}^{\dagger}\\\
n_{R}^{c\dagger}\\\ n_{D_{L}}^{\dagger}\\\
n_{D_{R}}^{c\dagger}\end{array}\right)=\mathcal{M}\left(\begin{array}[]{c}n_{L}\\\
n_{R}^{c}\\\ n_{D_{L}}\\\ n_{D_{R}}^{c}\end{array}\right),$ (16)

where

$\mathcal{M}=\left(\begin{array}[]{cccc}0&\tilde{m}_{N}&0&g_{R}^{*}\\\
\tilde{m}_{N}&0&g_{L}&0\\\ 0&g_{L}&M_{L}&m_{D}\\\
g_{R}^{*}&0&m_{D}&M_{R}\end{array}\right).$ (17)

Those can be transformed into the momentum space by the Fourier transform:

$\left(\begin{array}[]{c}n_{L}(t,\vec{x})\\\ n_{R}^{c}(t,\vec{x})\\\
n_{D_{L}}(t,\vec{x})\\\
n_{D_{R}}^{c}(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)\,,$ (18)

where

$\left(e_{\vec{k}}^{s}\right)_{1}=\sqrt{\frac{1}{2}\left(1+\frac{sk^{3}}{|\vec{k}|}\right)},\quad\left(e_{\vec{k}}^{s}\right)_{2}=se^{i\theta_{\vec{k}}}\sqrt{\frac{1}{2}\left(1-\frac{sk^{3}}{|\vec{k}|}\right)}\,,$
(19)

and

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

Then we get

$\displaystyle\partial_{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}+g_{R}D_{R,-\vec{k}}^{s\dagger}\right)\,,$
(21) $\displaystyle\partial_{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}+g_{L}^{*}D_{L,-\vec{k}}^{s\dagger}\right)\,,$
$\displaystyle\partial_{t}D_{L,\vec{k}}^{s}$
$\displaystyle=iskD_{L,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(m_{D}^{*}D_{R,-\vec{k}}^{s\dagger}+M_{L}^{*}D_{L,-\vec{k}}^{s\dagger}+g_{L}^{*}n_{R,-\vec{k}}^{s\dagger}\right)\,,$
$\displaystyle\partial_{t}D_{R,\vec{k}}^{s}$
$\displaystyle=iskD_{R,\vec{k}}^{s}+ise^{-i\theta_{\vec{k}}}\left(m_{D}^{*}D_{L,-\vec{k}}^{s\dagger}+M_{R}^{*}D_{R,-\vec{k}}^{s\dagger}+g_{R}n_{L,-\vec{k}}^{s\dagger}\right)\,.$

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

###  A.2 Time evolution of two-point functions and the total/net baryon number

The total number densities of neutron can be read off from Eqs.(8) and (9 as:

$n_{n}+\bar{n}_{n}=\frac{1}{V}\int\frac{d^{3}k}{(2\pi)^{3}}\sum_{s=\pm}\left(1-\frac{s|\vec{k}|}{\omega_{k}}\left(\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+\langle
n_{R,\vec{k}}^{s\dagger}n_{R,\vec{k}}^{s}\rangle\right)+\frac{\tilde{m}_{N}}{\omega_{k}}\left(se^{i\theta_{\vec{k}}}\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+se^{-i\theta_{\vec{k}}}\langle
n_{L,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle\right)\right)$ (22)

where $V=\int d^{3}x$ is a volume of the system and

$\omega_{k}(t)=\sqrt{|\vec{k}|^{2}+\tilde{m}_{N}(t)^{2}}.$ (23)

The net baryon number density in Eq.(13) is rewritten in terms of the Fourier
transformed fields as

$\begin{aligned} n_{B}&=\frac{1}{V}\int d^{3}x\frac{1}{2}\left(\left\langle
n_{L}^{\dagger}\bar{\sigma}^{0}n_{L}\right\rangle-\left\langle
n_{L}\sigma^{0}n_{L}^{\dagger}\right\rangle-\left\langle
n_{R}^{c\dagger}\bar{\sigma}^{0}n_{R}^{c}\right\rangle+\left\langle
n_{R}^{c}\sigma^{0}n_{R}^{c\dagger}\right\rangle\right)\\\
&=\frac{1}{V}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2}\sum_{s=\pm}\left(\left\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\right\rangle-\left\langle
n_{L,-\vec{k}}^{s}n_{L,-\vec{k}}^{s\dagger}\right\rangle\right).\\\
&=\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}\,,$
(24)

The two-point functions relevant to this net baryon number density are
developed in time through the time evolution equations given in Eq.(21):

$\displaystyle\partial_{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_{R}^{\ast}\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (25)
$\displaystyle\partial_{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}^{\ast}\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g_{L}^{\ast}\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}^{\ast}\langle
n_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\langle
n_{R,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (26)
$\displaystyle\partial_{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}^{\ast}\langle
D_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\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_{D}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+M_{R}^{\ast}\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g_{R}\langle
n_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (27)
$\displaystyle\partial_{t}\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=ise^{-i\theta_{\vec{k}}}\left(\tilde{m}_{N}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\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_{D}\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+M_{L}\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g_{L}\langle
n_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (28)
$\displaystyle\partial_{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_{L}\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (29)
$\displaystyle\partial_{t}\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{D}\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+M_{R}\langle
D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+g_{R}\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}^{\ast}\langle
D_{R,-\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\langle
D_{R,-\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (30)
$\displaystyle\partial_{t}\langle
D_{R,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{D}\langle
D_{L,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle+M_{R}\langle
D_{R,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle+g_{R}^{\ast}\langle
n_{L,\vec{k}}^{s}D_{R,-\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (31)
$\displaystyle\partial_{t}\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(m_{D}\langle
D_{L,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+M_{R}\langle
D_{R,-\vec{k}}^{s}n_{L,\vec{k}}^{s}\rangle+g_{R}^{\ast}\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}^{\ast}\langle
D_{R,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\langle
D_{R,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (32)
$\displaystyle\partial_{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_{D}^{\ast}\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+M_{L}^{\ast}\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g_{L}^{\ast}\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}^{\ast}\langle
D_{L,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\langle
D_{L,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (33)
$\displaystyle\partial_{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}^{\ast}\langle
D_{R,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\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_{D}^{\ast}\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+M_{R}^{\ast}\langle
D_{R,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g_{R}\langle
n_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle\right)\,,$ (34)
$\displaystyle\partial_{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_{D}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+M_{R}^{\ast}\langle
D_{R,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+g_{R}\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_{D}^{\ast}\langle
D_{R,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle+M_{L}^{\ast}\langle
D_{R,-\vec{k}}^{s}D_{L,-\vec{k}}^{s\dagger}\rangle+g_{L}^{*}\langle
D_{R,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (35)
$\displaystyle\partial_{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_{R}\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_{D}^{\ast}\langle
D_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+M_{L}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g_{L}^{*}\langle
n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle\right)\,,$ (36)
$\displaystyle\partial_{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}^{\ast}\langle
n_{R,\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s}\rangle+g_{R}\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}^{\ast}\langle
n_{L,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\langle
n_{L,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (37)
$\displaystyle\partial_{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_{D}^{\ast}\langle
D_{R,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle+M_{R}^{\ast}\langle
D_{R,\vec{k}}^{s}D_{R,\vec{k}}^{s\dagger}\rangle+g_{R}\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_{D}^{\ast}\langle
D_{L,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle+M_{R}^{\ast}\langle
D_{R,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle+g_{R}\langle
n_{L,-\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s}\rangle\right)\,,$ (38)
$\displaystyle\partial_{t}\langle
D_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle$
$\displaystyle=ise^{i\theta_{\vec{k}}}\left(m_{D}\langle
D_{R,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+M_{L}\langle
D_{L,\vec{k}}^{s}n_{R,-\vec{k}}^{s}\rangle+g_{L}\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}^{\ast}\langle
D_{L,-\vec{k}}^{s\dagger}n_{L,\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\langle
D_{L,-\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (39)
$\displaystyle\partial_{t}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle$
$\displaystyle=-ise^{i\theta_{\vec{k}}}\left(m_{D}\langle
D_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle+M_{L}\langle
D_{L,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle+g_{L}\langle
n_{R,-\vec{k}}^{s}D_{L,\vec{k}}^{s}\rangle\right)+{\rm h.c.}\,,$ (40)
$\displaystyle\partial_{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}^{\ast}\langle
n_{L,-\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s}\rangle+g_{L}^{\ast}\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}^{\ast}\langle
n_{R,\vec{k}}^{s}n_{L,\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\langle
n_{R,\vec{k}}^{s}D_{L,\vec{k}}^{s\dagger}\rangle\right)\,,$ (41)
$\displaystyle\partial_{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_{D}^{\ast}\langle
D_{R,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+M_{L}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,\vec{k}}^{s}\rangle+g_{L}^{\ast}\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_{D}^{\ast}\langle
D_{L,-\vec{k}}^{s}D_{R,-\vec{k}}^{s\dagger}\rangle+M_{L}^{\ast}\langle
D_{L,-\vec{k}}^{s}D_{L,-\vec{k}}^{s\dagger}\rangle+g_{L}^{\ast}\langle
D_{L,-\vec{k}}^{s}n_{R,-\vec{k}}^{s\dagger}\rangle\right)\,,$ (42)
$\displaystyle\partial_{t}\langle
D_{L,\vec{k}}^{s\dagger}D_{R,\vec{k}}^{s}\rangle$
$\displaystyle=ise^{-i\theta_{\vec{k}}}\left(m_{D}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}D_{L,-\vec{k}}^{s\dagger}\rangle+M_{R}^{\ast}\langle
D_{L,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\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_{D}\langle
D_{R,-\vec{k}}^{s}D_{R,\vec{k}}^{s}\rangle+M_{L}\langle
D_{L,-\vec{k}}^{s}D_{R,\vec{k}}^{s}\rangle+g_{L}\langle
n_{R,-\vec{k}}^{s}D_{R,\vec{k}}^{s}\rangle\right)\,,$ (43)
$\displaystyle\partial_{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_{L}\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}^{\ast}\langle
n_{R,\vec{k}}^{s\dagger}n_{R,-\vec{k}}^{s\dagger}\rangle+g_{R}\langle
n_{R,\vec{k}}^{s\dagger}D_{R,-\vec{k}}^{s\dagger}\rangle\right)\,.$ (44)

###  A.3 Initial values

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 since the fields do not have any other
interactions except those due to $\langle\sigma\rangle$. It is convenient to
introduce a collective Majonara 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}\,.$ (45)

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}\,,$
(46)

their initial conditions at $t=t_{0}$ are then set as

$\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}$ (47)

where $\mathcal{M}_{d}$ is a diagonalized mass matrix related to the original
mass matrix $\mathcal{M}$ in Eq.(17) through a unitary transformation by $U$
as

$\mathcal{M}_{d}\equiv U^{T}\mathcal{M}U\,.$ (48)

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

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

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

$\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}$
(50)

## References

  * (1) S. Ipek and T. M. P. Tait, Phys. Rev. Lett. 122 (2019) no.11, 112001 doi:10.1103/PhysRevLett.122.112001 [arXiv:1811.00559 [hep-ph]].
  * (2) G. Servant, Phys. Rev. Lett. 113 (2014) no.17, 171803 doi:10.1103/PhysRevLett.113.171803 [arXiv:1407.0030 [hep-ph]].
  * (3) D. Croon, J. N. Howard, S. Ipek and T. M. P. Tait, Phys. Rev. D 101 (2020) no.5, 055042 doi:10.1103/PhysRevD.101.055042 [arXiv:1911.01432 [hep-ph]].
  * (4) A. D. Dolgov and D. P. Kirilova, Sov. J. Nucl. Phys. 51 (1990), 172-177 JINR-E2-89-321.
  * (5) J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 42 (1990), 2491-2504 doi:10.1103/PhysRevD.42.2491
  * (6) L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. Lett. 73 (1994), 3195-3198 doi:10.1103/PhysRevLett.73.3195 [arXiv:hep-th/9405187 [hep-th]].
  * (7) Y. Shtanov, J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 51 (1995), 5438-5455 doi:10.1103/PhysRevD.51.5438 [arXiv:hep-ph/9407247 [hep-ph]].
  * (8) L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56 (1997), 3258-3295 doi:10.1103/PhysRevD.56.3258 [arXiv:hep-ph/9704452 [hep-ph]].
  * (9) M. A. Amin, M. P. Hertzberg, D. I. Kaiser and J. Karouby, Int. J. Mod. Phys. D 24 (2014), 1530003 doi:10.1142/S0218271815300037 [arXiv:1410.3808 [hep-ph]].
  * (10) K. D. Lozanov, [arXiv:1907.04402 [astro-ph.CO]].
  * (11) S. Enomoto, C. Cai, Z. H. Yu, H. H. Zhang, “Matter-antimatter asymmetry in preheating,” AAPPS Bull. 30 (2020) no. 5, 45–48. http://aappsbulletin.org/myboard/read.php?id=125&Page=2&Board=focus&FindIt=&FindText=
  * (12) K. Funakubo, A. Kakuto, S. Otsuki and F. Toyoda, Prog. Theor. Phys. 105 (2001), 773-788 doi:10.1143/PTP.105.773 [arXiv:hep-ph/0010266 [hep-ph]].
  * (13) R. Rangarajan and D. V. Nanopoulos, Phys. Rev. D 64 (2001), 063511 doi:10.1103/PhysRevD.64.063511 [arXiv:hep-ph/0103348 [hep-ph]].
  * (14) L. Pearce, L. Yang, A. Kusenko and M. Peloso, Phys. Rev. D 92 (2015) no.2, 023509 doi:10.1103/PhysRevD.92.023509 [arXiv:1505.02461 [hep-ph]].
  * (15) S. Enomoto and T. Matsuda, Int. J. Mod. Phys. A 33 (2018) no.25, 1850146 doi:10.1142/S0217751X18501464 [arXiv:1707.05310 [hep-ph]].
  * (16) Y. P. Wu, L. Yang and A. Kusenko, JHEP 12 (2019), 088 doi:10.1007/JHEP12(2019)088 [arXiv:1905.10537 [hep-ph]].
  * (17) S. Enomoto, C. Cai, Z. H. Yu and H. H. Zhang, Eur. Phys. J. C 80 (2020) no.12, 1098 doi:10.1140/epjc/s10052-020-08672-7 [arXiv:2005.08037 [hep-ph]].
  * (18) S. M. Lee, K. y. Oda and S. C. Park, JHEP 03 (2021), 083 doi:10.1007/JHEP03(2021)083 [arXiv:2010.07563 [hep-ph]].
  * (19) E. W. Kolb, A. D. Linde and A. Riotto, Phys. Rev. Lett. 77 (1996), 4290-4293 doi:10.1103/PhysRevLett.77.4290 [arXiv:hep-ph/9606260 [hep-ph]].
  * (20) E. W. Kolb, A. Riotto and I. I. Tkachev, Phys. Lett. B 423 (1998), 348-354 doi:10.1016/S0370-2693(98)00134-8 [arXiv:hep-ph/9801306 [hep-ph]].
  * (21) G. F. Giudice, M. Peloso, A. Riotto and I. Tkachev, JHEP 08 (1999), 014 doi:10.1088/1126-6708/1999/08/014 [arXiv:hep-ph/9905242 [hep-ph]].
  * (22) M. Peloso and L. Sorbo, JHEP 05 (2000), 016 doi:10.1088/1126-6708/2000/05/016 [arXiv:hep-ph/0003045 [hep-ph]].
  * (23) J. Garcia-Bellido, D. G. Figueroa and J. Rubio, Phys. Rev. D 79 (2009), 063531 doi:10.1103/PhysRevD.79.063531 [arXiv:0812.4624 [hep-ph]].
  * (24) P. Adshead and E. I. Sfakianakis, JCAP 11 (2015), 021 doi:10.1088/1475-7516/2015/11/021 [arXiv:1508.00891 [hep-ph]].
  * (25) Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, JCAP 02 (2017), 045 doi:10.1088/1475-7516/2017/02/045 [arXiv:1609.05209 [hep-ph]].
  * (26) P. Adshead, L. Pearce, M. Peloso, M. A. Roberts and L. Sorbo, JCAP 06 (2018), 020 doi:10.1088/1475-7516/2018/06/020 [arXiv:1803.04501 [astro-ph.CO]].
  * (27) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 06 (2009), 088 doi:10.1088/1126-6708/2009/06/088 [arXiv:0903.4155 [hep-lat]].
  * (28) S. Borsanyi et al. [Wuppertal-Budapest], J. Phys. Conf. Ser. 316 (2011), 012020 doi:10.1088/1742-6596/316/1/012020 [arXiv:1109.5032 [hep-lat]].
  * (29) H. T. Ding, F. Karsch and S. Mukherjee, Int. J. Mod. Phys. E 24 (2015) no.10, 1530007 doi:10.1142/S0218301315300076 [arXiv:1504.05274 [hep-lat]].
  * (30) A. Bazavov et al. [HotQCD], Phys. Lett. B 795 (2019), 15-21 doi:10.1016/j.physletb.2019.05.013 [arXiv:1812.08235 [hep-lat]].
  * (31) H. T. Ding, Nucl. Phys. A 1005 (2021), 121940 doi:10.1016/j.nuclphysa.2020.121940 [arXiv:2002.11957 [hep-lat]].
  * (32) A. Bazavov, N. Brambilla, H. T. Ding, P. Petreczky, H. P. Schadler, A. Vairo and J. H. Weber, Phys. Rev. D 93 (2016) no.11, 114502 doi:10.1103/PhysRevD.93.114502 [arXiv:1603.06637 [hep-lat]].
  * (33) H. T. Ding, PoS LATTICE2016 (2017), 022 doi:10.22323/1.256.0022 [arXiv:1702.00151 [hep-lat]].
  * (34) R. D. Pisarski and F. Wilczek, Phys. Rev. D 29 (1984), 338-341 doi:10.1103/PhysRevD.29.338
  * (35) S. Iso, P. D. Serpico and K. Shimada, Phys. Rev. Lett. 119 (2017) no.14, 141301 doi:10.1103/PhysRevLett.119.141301 [arXiv:1704.04955 [hep-ph]].
  * (36) T. Hambye, A. Strumia and D. Teresi, JHEP 08 (2018), 188 doi:10.1007/JHEP08(2018)188 [arXiv:1805.01473 [hep-ph]].
  * (37) L. Sagunski, P. Schicho and D. Schmitt, [arXiv:2303.02450 [hep-ph]].
  * (38) J. E. Kim, S. J. Kim and S. Nam, [arXiv:1803.03517 [hep-ph]].
  * (39) J. E. Kim and S. J. Kim, Phys. Lett. B 783 (2018), 357-363 doi:10.1016/j.physletb.2018.07.020 [arXiv:1804.05173 [hep-ph]].
  * (40) J. E. Kim and S. J. Kim, PoS CORFU2018 (2019), 052 doi:10.22323/1.347.0052 [arXiv:1902.02766 [hep-ph]].
  * (41) P. A. Zyla et al. [Particle Data Group], PTEP 2020 (2020) no.8, 083C01 doi:10.1093/ptep/ptaa104
  * (42) A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5 (1967), 32-35 doi:10.1070/PU1991v034n05ABEH002497
  * (43) S. Enomoto and T. Matsuda, Phys. Rev. D 99 (2019) no.3, 036005 doi:10.1103/PhysRevD.99.036005 [arXiv:1811.06197 [hep-th]].
  * (44) T. Bringmann, J. M. Cline and J. M. Cornell, Phys. Rev. D 99 (2019) no.3, 035024 doi:10.1103/PhysRevD.99.035024 [arXiv:1810.08215 [hep-ph]].