# A New Mechanism for Generating Particle Number Asymmetry through Interactions

Authors: Takuya Morozumi; Keiko I. Nagao; Apriadi Salim Adam; Hiroyuki Takata
Journal: Advances in High Energy Physics
 (2019)

---

## Abstract

A new mechanism for generating particle number asymmetry (PNA) has been developed. This mechanism is realized with a Lagrangian including a complex scalar field and a neutral scalar field. The complex scalar carriesU(1) charge which is associated with the PNA. It is written in terms of the condensation and Green’s function, which is obtained with two-particle irreducible (2PI) closed time path (CTP) effective action (EA). In the spatially flat universe with a time-dependent scale factor, the time evolution of the PNA is computed. We start with an initial condition where only the condensation of the neutral scalar is nonzero. The initial condition for the fields is specified by a density operator parameterized by the temperature of the universe. With the above initial conditions, the PNA vanishes at the initial time and later it is generated through the interaction between the complex scalar and the condensation of the neutral scalar. We investigate the case that both the interaction and the expansion rate of the universe are small and include their effects up to the first order of the perturbation. The expanding universe causes the effects of the dilution of the PNA, freezing interaction, and the redshift of the particle energy. As for the time dependence of the PNA, we found that PNA oscillates at the early time and it begins to dump at the later time. The period and the amplitude of the oscillation depend on the mass spectrum of the model, the temperature, and the expansion rate of the universe.


---

## Body

## 1. Introduction

The origin of BAU has long been a question of great interest in explaining why there is more baryon than antibaryon in nature. Big bang nucleosynthesis (BBN) [1] and cosmic microwave background [2] measurements give the BAU as η≡nB/s≃10-10, where nB is the baryon number density and s is the entropy density. In order to address this issue, many different models and mechanisms have been proposed [3–7]. The mechanisms discussed in the literature satisfy the three Sakharov conditions [3], namely, (i) baryon number (B) violation, (ii) charge (C) and charge-parity (CP) violations, and (iii) a departure from the thermal equilibrium. For reviews of different types of models and mechanisms, see, for example, [8–10]. Recently, the variety of the method for the calculation of BAU has been also developed [11–13].In the present paper, we further extend the model of scalar fields [14] so that it generates the PNA through interactions. In many of the previous works, the mechanism generating BAU relies on the heavy particle decays. Another mechanism uses U(1) phase of the complex scalar field [6]. In this work, we develop a new mechanism to generate PNA. The new feature of our approach is briefly explained below.The model which we have proposed [15] consists of a complex scalar field and a neutral scalar field. The PNA is related to the U(1) current of the complex field. In our model, the neutral scalar field has a time-dependent expectation value which is called condensation. In the new mechanism, the oscillating condensation of the neutral scalar interacts with the complex scalar field. Since the complex scalar field carries U(1) charge, the interactions with the condensation of the neutral scalar generate PNA. The interactions break U(1) symmetry as well as charge conjugation symmetry. At the initial time, the condensation of the neutral scalar is nonzero. We propose a way which realizes such initial condition.As for the computation of the PNA, we use 2PI formalism combined with density operator formulation of quantum field theory [16]. The initial conditions of the quantum fields are specified with the density operator. The density operator is parameterized by the temperature of the universe at the initial time. We also include the effect of the expansion of the universe. It is treated perturbatively and the leading order term which is proportional to the Hubble parameter at the initial time is considered. With this method, the time dependence of the PNA is computed and the numerical analysis is carried out. The dependence, especially, on the various parameters of the model such as masses and strength of interactions is investigated. We also study the dependence on the temperature and the Hubble parameter at the initial time. We first carry out the numerical simulation without specifying the unit of parameter sets. Later, in a radiation dominated era, we specify the unit of the parameters and estimate the numerical value of the PNA over entropy density.This paper is organized as follows. In Section2, we introduce our model with CP and particle number violating interactions. We also specify the density operator as the initial state. In Section 3, we derive the equation of motion for Green’s function and field by using 2PI CTP EA formalism. We also provide the initial condition for Green’s function and field. In Section 4, using the solution of Green’s function and field, we compute the expectation value of the PNA. Section 5 provides the numerical study of the time dependence of the PNA. We will also discuss the dependence on the parameters of the model. Section 6 is devoted to conclusion and discussion. In Appendix A, we introduce a differential equation which is a prototype for Green’s function and field equations. Applying the solutions of the prototype, we obtain the solutions for both Green’s function and field equations. In Appendices B–D, the useful formulas to obtain the PNA for nonvanishing Hubble parameter case are derived.


## 2. A Model with CP and Particle Number Violating Interaction

In this section, we present a model which consists of scalar fields [15]. It has both CP and particle number violating features. As an initial statistical state for scalar fields, we employ the density operator for thermal equilibrium.Let us start by introducing a model consisting of a neutral scalar,N, and a complex scalar, ϕ. The action is given by(1)S=∫d4x-gLfree+Lint,Lfree=gμν∇μϕ†∇νϕ-mϕ2ϕ2+12∇μN∇μN-MN22N2+B22ϕ2+ϕ†2+α22ϕ2+h.c.R+α3ϕ2R,Lint=Aϕ2N+A∗ϕ†2N+A0ϕ2N,where gμν is the metric and R is the Riemann curvature. With this Lagrangian, we aim to produce the PNA through the soft-breaking terms of U(1) symmetry whose coefficients are denoted by A and B2. One may add the quartic terms to the Lagrangian which are invariant under the U(1) symmetry. Though those terms preserve the stability of the potential for large field configuration and are also important for the renormalizability, we assume they do not lead to the leading contribution for the generation of the PNA. We also set the coefficients of the odd power terms for Nn(n=1,3) zero in order to obtain a simple oscillating behavior for the time dependence of the condensation of N. We assume that our universe is homogeneous for space and employ the Friedmann-Lemaître-Robertson-Walker metric,(2)gμν=1,-a2x0,-a2x0,-a2x0,where a(x0) is the scale factor at time x0. Correspondingly the Riemann curvature is given by(3)Rx0=6a¨x0ax0+a˙x0ax02.In (1), the terms proportional to A, B, and α2 are the particle number violating interactions. In general, only one of the phases of those parameters can be rotated away. Throughout this paper, we study the special case that B and α2 are real numbers and A is a complex number. Since only A is a complex number, it is a unique source of the CP violation.We rewrite all the fields in terms of real scalar fields,ϕi (i=1,2,3), defined as(4)ϕ=ϕ1+iϕ22,N=ϕ3.With these definitions, the free part of the Lagrangian is rewritten as(5)Lfree=12-ggμν∇μϕi∇νϕi-m~i2x0ϕi2,where the kinetic term is given by(6)gμν∇μϕi∇νϕi=∂ϕi∂x0∂ϕi∂x0-1ax02∂ϕi∂xj∂ϕi∂xj,and their effective masses are given as follows:(7)m~12x0=mϕ2-B2-α2+α3Rx0,(8)m~22x0=mϕ2+B2+α2-α3Rx0,(9)m~32x0=mN2.Nonzero B2 or α2 leads to the nondegenerate mass spectrum for ϕ1 and ϕ2. The interaction Lagrangian is rewritten with a totally symmetric coefficient Aijk,(10)Lint=∑ijk=1313Aijkϕiϕjϕkwith i,j,k=1,2,3. The nonzero components of Aijk are written with the couplings for cubic interaction, A and A0, as shown in Table 1. We also summarize the qubic interactions and their properties according to U(1) symmetry and CP symmetry.Table 1
The cubic interactions and their properties.
Cubic interaction coupling Property A113=A02+Re.(A) –  A223=A02-Re.(A) –  A113-A223=2Re.(A) U ( 1 ) violation  A123=-Im.(A) U ( 1 ), CP violationNöether current related to theU(1) transformation is(11)jμx=i2ϕ†∂μ⃡ϕ-ϕ∂μ⃡ϕ†.In terms of real scalar fields, the Nöether current alters into(12)jμ=12ϕ2∂μ⃡ϕ1-ϕ1∂μ⃡ϕ2.The ordering of the operators in (11) is arranged so that it is Hermite and the particle number operator,(13)Qx0=∫d3x-gj0x,has a normal ordered expression. Then, in the vanishing limit of interaction terms and particle number violating terms, the vacuum expectation value of the particle number vanishes. With the above definition, j0(x) is the PNA per unit comoving volume. The expectation value of the PNA is written with a density operator,(14)j0x=Trj0xρt0.Note that, the PNA is a Heisenberg operator and ρ(t0) is a density operator which specifies the state at the initial time x0=t0. In this work, we use the density operator with zero chemical potential. It is specifically given by(15)ρt0=e-βHTre-βH,where β denotes inverse temperature, 1/T, and H is a Hamiltonian which includes linear term of fields,(16)H=12∑i=13∫d3xat03πϕiπϕi+∇ϕi·∇ϕiat02+m~i2ϕi-vi2,where vi is a constant. The linear term of fields in (16) is prepared for the nonzero expectation value of fields. Note that the density operator in (15) is not exactly the same as the thermal equilibrium one since, in the Hamiltonian, the interaction parts are not included. Since we assume three-dimensional space is translational invariant, then the expectation value of the PNA depends on time x0 and the initial time t0. As we will show later, the nonzero expectation value for the field ϕ3 leads to the time-dependent condensation which is the origin of the nonequilibrium time evolution of the system.Below we consider the matrix element of the density operator given in (15). We start with the following density operator for one component real scalar field as an example,(17)ρt0=e-βHexampleTre-βHexample,(18)Hexample=12∫d3xat03πϕπϕ+∇ϕ·∇ϕat02+m~2ϕ-v2.The above Hamiltonian is obtained from that of (16) by keeping only one of the real scalar fields. The matrix element of the initial density operator in (17) is written in terms of the path integral form of the imaginary time formalism given as(19)ϕ1ρt0ϕ2=∫ϕ0=ϕ2,ϕβ=ϕ1dϕe-SEexampleϕ∫dϕ1∫ϕ0=ϕ1,ϕβ=ϕ1dϕe-SEexampleϕ,where SEexample is an Euclidean action which corresponds to the Hamiltonian in (18) and it is given by(20)SEexampleϕx,u=12∫0βdu∫d3x∂ϕ∂u2+∇ϕ·∇ϕat02+m~2ϕ-v2.After carrying out the path integral, the density matrix is written with SEexample[ϕcl(x,u)] which is the action for the classical orbit ϕcl satisfying the boundary conditions, ϕcl(u=0)=ϕ2,ϕcl(u=β)=ϕ1. It is given as the functional of the boundary fields ϕi(i=1,2) and vacuum expectation value v as(21)ϕ1ρt0ϕ2=exp⁡-SEclexampleϕ1,ϕ2∫dϕ1exp⁡-SEclexampleϕ1,ϕ1,where SEclexample[ϕ1,ϕ2]≃SEexample[ϕcl(x,u)] is given by(22)SEclexampleϕ1,ϕ2=-at062∫d3k2π3∑i=12ϕikκii-kϕi-k-ϕ1kκ12-kϕ2-k-ϕ2kκ21-kϕ1-k+at06∑i=12ϕi0κii0-ϕ10κ120-ϕ20κ210v.In the above expression, we drop the terms which are proportional to v2 because they do not contribute to the normalized density matrix. κbd(k) is defined as [14](23)κ11k=κ22k≔-1at03ωkcosh⁡βωksinh⁡βωk,κ12k=κ21k≔-1at03ωksinh⁡βωk,ωk≔k2at02+m~2.Using the above definitions, one can write the density matrix in (21) as the following form,(24)ϕ1ρt0ϕ2=Nexp⁡∫-gxd4xJaxcabϕbx+12∫d4xd4y-gxϕaxcabKbdx,ycdeϕey-gy,with ϕ1(t0)=ϕ1 and ϕ2(t0)=ϕ2. The upper indices a and b are 1 or 2. cab is the metric of CTP formalism [16] and c11=-c22=1 and c12=c21=0. In the above expression, the source terms J and K do not vanish only at the initial time t0 and they are given by(25)Jbx≔-iδx0-t0jb,jb=-at03κbdk=0cdeve,(26)Kbdx,y≔-iδx0-t0δy0-t0κbdx-y,κabx=∫d3k2π3κabke-ik·x,where κab(k) is given in (23). In (24), N is a normalization constant which is given as(27)1N=∫dϕ1∫dϕ2δϕ1-ϕ2exp⁡∫-gxd4xJaxcabϕbx+12∫d4xd4y-gxϕaxcabKbdx,ycdeϕey-gy.


## 3. Two-Particle Irreducible Closed Time Path Effective Action

In this section, we derive the equations of motion, i.e., the Schwinger-Dyson equations (SDEs) for both Green’s function and field. SDEs are obtained by taking the variation of 2PI EA with respect to fields and Green’s functions, respectively. In addition, we also provide the initial condition for Green’s function and field to solve SDEs.


### 3.1. 2PI Formalism in Curved Space-Time

2PI CTP EA in curved space-time has been investigated in [17] and their formulations can be applied to the present model. In 2PI formalism, one introduces nonlocal source term denoted as K and local source term denoted as J,(28)eiWJ,K=∫dϕexp⁡iSϕ,g+∫-gxd4xJiaxcabϕibx+12∫d4xd4y-gx×ϕiaxcabKijbdx,ycdeϕjey-gy,where i,j=1 or 2 and S[ϕ,g] is given by(29)Sϕ,g=∫d4x-gx12cabgμν∇μϕia∇νϕib-m~ii2ϕiaϕib+13DabcAijkϕiaϕjbϕkc,where D111=-D222=1 and the other components are zero. The upper indices of the field and the source terms distinguish two different time paths in closed time path formalism [16]. One can define the mean fields ϕ¯ia and Green’s function by taking the functional derivative with respect to the source terms J and K, respectively,(30)ϕ¯iax=cab-gxδWJ,KδJibx,(31)ϕ¯iaxϕ¯jey+Gijaex,y=2cab-gxδWJ,KδKijbdx,ycde-gy.If one sets the source terms to be the ones given in (25) and (26), one can show that the expectation value of the product of the field operators with the initial density operator is related to the Green function and mean fields. Definitely, we can prove the following relations,(32)ϕ¯iax=∫dϕ1∫dϕ2ϕ2ϕixϕ1exp⁡-SEclϕ1,ϕ2∫dϕ1∫dϕ2δϕ1-ϕ2exp⁡-SEclϕ1,ϕ2(33)=Trϕixρt0,(34)Gij12x,y=∫∫dϕ1dϕ2ϕ2ΦjyΦixϕ1e-SEclϕ1,ϕ2∫∫dϕ1dϕ2δϕ2-ϕ1e-SEclϕ1,ϕ2(35)=Trϕjyϕixρt0-ϕ¯j2yϕ¯i1x,with ϕ¯a(x)=ϕ¯(x) and Φ(x) is a Heisenberg operator which has form as Φ(x0,x)≡ϕ(x0,x)-ϕ¯(x0). With (35), one can write the expectation value of the current as the sum of the contribution from Green’s function and the current of the mean fields. Then (14) alters into(36)j0x=Re∂∂x0-∂∂y0G1212x,yy→x+ϕ¯22x∂0⃡ϕ¯11x,where we have used (12) and the following relations(37)Gijab∗x,y=τ1acGijcdx,yτ1db,(38)ϕ¯ia∗x=τ1abϕ¯ibx,where τ is the Pauli matrix.The Green’s functions and expectation value of fields are derived as solutions of the SDEs which are obtained with 2PI EA. The 2PI EA is related to the generating functionalW[J,K] by Legendre transformation as [18, 19](39)Γ2G,ϕ¯,g=WJ,K-∫d4x-gxJiaxcabϕ¯ibx-12∫d4x∫d4y-gxcabKijbdx,ycdeϕ¯iaxϕ¯jey+Gijaex,y-gy.Let us write the 2PI EAΓ2 in our model, in which we only keep the interaction term up to the first order of cubic interaction, Aijk. It is given as(40)Γ2G,ϕ¯,g=Sϕ¯,g+12∫d4x∫d4yδ2Sϕ¯,gδϕ¯iaxδϕ¯jbyGijabx,y+i2TrLnG-1,where S[ϕ¯,g] is the action written in terms of mean fields as(41)Sϕ¯,g=∫d4x-gx-ϕ¯ia12cab□+m~ii2ϕ¯ib+13DabcAijkϕ¯iaϕ¯jbϕ¯kc+12∫d4x-gxδx0-T-δx0-t0ϕ¯iacabϕ¯˙ib.In (40), the interactions are included in the first term as well as in the second term. In the action above, we have also taken into account the surface term at the boundary which corresponds to the last term of (41). T and t0 in (41) are the upper bound and the lower bound of the time integration, respectively.


### 3.2. Schwinger Dyson Equations

Now let us derive SDEs for both Green’s function and field. These equations can be obtained by taking the variation of the 2PI EA,Γ2, with respect to the scalar field ϕ¯ and Green’s function G.In the following, we first derive SDEs for the field. The variation of the 2PI EA in (39) with respect to the scalar field ϕ¯ leads to(42)1-gxδΓ2δϕ¯iax=-cabJibx-∫d4zcabKijbcx,zccd-gzϕ¯jdz.Using (25) and (26), one computes the right-hand side of the above equation as(43)cabJibx+∫d4zcabKijbcx,zccd-gzϕ¯jdz=-iδx0-t0vim~itanh⁡βm~i2+cabκiibck=0ccdat03vid=0,where we have used κab(k) given in (23). The left-hand side of Eq. (42) is computed using (40) and one obtains the following equation of motion of the scalar field ϕ¯,(44)δij□+m~ij2ϕ¯jdx=cdaDabcAijkϕ¯jbxϕ¯kcx+Gjkbcx,x,where the Laplacian of Friedman-Lemaître-Robertson-Walker metric is given by(45)□=∇μ∇μ=∂2∂x02-1ax02∇·∇+3a˙a∂∂x0.Next, the equation of motion for Green’s function is derived in the following way. The variation of the 2PI EA in (39) with respect to Green’s function G leads to(46)δΓ2δGijabx,y=-12cac-gxKijcdx,y-gycdb.The left-hand side of the above equation is obtained by taking variation of Eq. (40) with respect to Green’s function as,(47)δΓ2δGijabx,y=-i2G-1jibay,x+12δ2Sϕ¯,gδϕ¯iaxδϕ¯jby,where the second term of above expression is computed using action in (41). Taking all together (46) and (47), one obtains the following two differential equations for Green’s function,(48)□→x+m~i2Gijabx,y=-iδijcab-gxδx-y+2cadDdceAiklϕ¯l,xeGkj,xycb+∫d4zKikaex,z-gzcefGkjfbz,y,(49)□→y+m~j2Gijabx,y=-iδijδx-ycab-gy+2Gik,xyacDcefAkjlϕ¯l,yfceb+∫d4zGikaex,zcef-gzKkjfbz,y,where □x=∇xμ∇μx and □y=∇yμ∇μy.Next, we rescale Green’s function, field and coupling constant of interaction as follows:(50)ϕ¯x0≕at0ax03/2φ^x0,(51)Gx0,y0,k≕at0ax03/2G^x0,y0,kat0ay03/2,(52)A^x0≔at0ax03/2A,where at0 stands for the initial value for the scale factor and we have defined at0≔a(t0) and we have used Fourier transformation for Green’s function as(53)Gx0,y0,k=∫d3rGx0,r,y0,0eik·r.By using these new definitions, SDEs in (44) are written as(54)∂2∂x02+m¯i2x0φ^idx0=cdaDabcA^ijkx0φ^jbx0φ^kcx0+G^jkbcx,x.Next SDEs for the rescaled Green’s function in (48) and (49) are written as(55)∂2∂x02+Ωi,k2x0G^ij,x0y0abk=2cadDdceA^ikl,x0φ^l,x0eG^kj,x0y0cbk-iδijδx0y0cabat03-iδt0x0κikaekat03cefG^kj,t0y0fbk,(56)∂2∂y02+Ωi,k2y0G^ij,x0y0abk=2G^ik,x0y0ackDcefA^kjl,y0φ^l,y0fceb-iδijδx0y0cabat03-iG^ik,x0t0aekcefat03κkjfbkδt0y0,where we have defined(57)Ωi,k2x0≔k2ax02+m¯i2x0,(58)m¯i2x0≔m~i2x0-32a¨x0ax0-34a˙x0ax02.Note that the first derivative with respect to time which is originally presented in the expression of Laplacian (45) is now absent in the expression of SDEs for the rescaled fields and Green’s functions.


### 3.3. The Initial Condition for Green’s Function and Field

In this subsection, the initial conditions for Green’s function and field are determined. For simplicity, let us look back to example model for one real scalar field. We first compute the initial condensation of the fieldϕ¯(t0)=φ^(t0) (see (50)). Using (32) and setting x0=t0, we compute it as follows:(59)ϕ¯t0≡ϕt0,x≔∫dϕϕt0,xexp⁡-SEclexampleϕ,ϕ∫dϕexp⁡-SEclexampleϕ,ϕ=∫dϕϕxexp⁡-1/2∫d3xd3yϕxDx-yϕy+2at03∫d3xϕxj1∫dϕexp⁡-1/2∫d3xd3yϕxDx-yϕy+2at03∫d3xϕxj1,where we have computed the last term of (22) using (25) and D(r) is defined as [14](60)Dr=2at03∫d3k2π3ωkcosh⁡βωk-1sinh⁡βωke-ir·k.To proceed the calculation, we denote J(x) as J(x)=2at03j1. Then the initial condensation of field 〈ϕ(t0,x)〉 is given by(61)ϕt0,x=∫dϕ′xϕ′x+∫d3zD-1x-zJzexp-1/2∫d3xd3yϕ′xDx-yϕ′y∫dϕ′xexp-1/2∫d3xd3yϕ′xDx-yϕ′y=v,where we have defined ϕ′(x)=ϕ(x)-∫d3zD-1(z-x)J(z) and D-1(z-x) satisfies(62)∫d3xD-1z-xDx-y=δz-y.Next we will compute the initial condition for Green’s functionG(t0,t0,k)=G^(t0,t0,k) (see (51)). By using (34) and setting x0=t0 and y0=t0, one computes it as follows:(63)Gabt0,x,t0,y=∫∫dϕ1dϕ2ϕ2Φt0,yΦt0,xϕ1exp⁡-SEclexampleϕ1,ϕ2∫dϕexp⁡-SEclexampleϕ,ϕ=∫dϕ′ϕ′yϕ′xexp⁡-1/2∫d3xd3yϕ′xDx-yϕ′y∫dϕ′exp⁡-1/2∫d3xd3yϕ′xDx-yϕ′y,=D-1y-x.Using (53) and (60), (63) becomes(64)Gabt0,t0,k=D-1k=12ωkat03sinh⁡βωkcosh⁡βωk-1.The above results with the example model can be extended to our model and we summarize them as(65)G^ij,t0t0abk=δij12ωikat03sinh⁡βωikcosh⁡βωik-1,(66)φ^iat0=vi.Next we derive the time derivative of the field and Green’s function at the initial timet0. First we integrate the field equation in (54) with respect to time. By setting x0=t0, we obtain(67)∂φ^i,x0∂x0x0=t0=0.Similarly, we integrate (55) with respect to time x0. By setting both x0 and y0 equal to t0, we obtain the following initial condition(68)limx0→t0⁡∂∂x0G^ij,x0t0abk=-iδijcabat03-iκikaekat03cefG^kj,t0t0fbk.Finally, we integrate (56) with respect to time y0. By setting both x0 and y0 equal to t0, we obtain another initial condition,(69)limy0→t0⁡∂∂y0G^ij,t0y0abk=-iδijcabat03-iG^ik,t0t0aekcefat03κkjfbk.


## 3.1. 2PI Formalism in Curved Space-Time

2PI CTP EA in curved space-time has been investigated in [17] and their formulations can be applied to the present model. In 2PI formalism, one introduces nonlocal source term denoted as K and local source term denoted as J,(28)eiWJ,K=∫dϕexp⁡iSϕ,g+∫-gxd4xJiaxcabϕibx+12∫d4xd4y-gx×ϕiaxcabKijbdx,ycdeϕjey-gy,where i,j=1 or 2 and S[ϕ,g] is given by(29)Sϕ,g=∫d4x-gx12cabgμν∇μϕia∇νϕib-m~ii2ϕiaϕib+13DabcAijkϕiaϕjbϕkc,where D111=-D222=1 and the other components are zero. The upper indices of the field and the source terms distinguish two different time paths in closed time path formalism [16]. One can define the mean fields ϕ¯ia and Green’s function by taking the functional derivative with respect to the source terms J and K, respectively,(30)ϕ¯iax=cab-gxδWJ,KδJibx,(31)ϕ¯iaxϕ¯jey+Gijaex,y=2cab-gxδWJ,KδKijbdx,ycde-gy.If one sets the source terms to be the ones given in (25) and (26), one can show that the expectation value of the product of the field operators with the initial density operator is related to the Green function and mean fields. Definitely, we can prove the following relations,(32)ϕ¯iax=∫dϕ1∫dϕ2ϕ2ϕixϕ1exp⁡-SEclϕ1,ϕ2∫dϕ1∫dϕ2δϕ1-ϕ2exp⁡-SEclϕ1,ϕ2(33)=Trϕixρt0,(34)Gij12x,y=∫∫dϕ1dϕ2ϕ2ΦjyΦixϕ1e-SEclϕ1,ϕ2∫∫dϕ1dϕ2δϕ2-ϕ1e-SEclϕ1,ϕ2(35)=Trϕjyϕixρt0-ϕ¯j2yϕ¯i1x,with ϕ¯a(x)=ϕ¯(x) and Φ(x) is a Heisenberg operator which has form as Φ(x0,x)≡ϕ(x0,x)-ϕ¯(x0). With (35), one can write the expectation value of the current as the sum of the contribution from Green’s function and the current of the mean fields. Then (14) alters into(36)j0x=Re∂∂x0-∂∂y0G1212x,yy→x+ϕ¯22x∂0⃡ϕ¯11x,where we have used (12) and the following relations(37)Gijab∗x,y=τ1acGijcdx,yτ1db,(38)ϕ¯ia∗x=τ1abϕ¯ibx,where τ is the Pauli matrix.The Green’s functions and expectation value of fields are derived as solutions of the SDEs which are obtained with 2PI EA. The 2PI EA is related to the generating functionalW[J,K] by Legendre transformation as [18, 19](39)Γ2G,ϕ¯,g=WJ,K-∫d4x-gxJiaxcabϕ¯ibx-12∫d4x∫d4y-gxcabKijbdx,ycdeϕ¯iaxϕ¯jey+Gijaex,y-gy.Let us write the 2PI EAΓ2 in our model, in which we only keep the interaction term up to the first order of cubic interaction, Aijk. It is given as(40)Γ2G,ϕ¯,g=Sϕ¯,g+12∫d4x∫d4yδ2Sϕ¯,gδϕ¯iaxδϕ¯jbyGijabx,y+i2TrLnG-1,where S[ϕ¯,g] is the action written in terms of mean fields as(41)Sϕ¯,g=∫d4x-gx-ϕ¯ia12cab□+m~ii2ϕ¯ib+13DabcAijkϕ¯iaϕ¯jbϕ¯kc+12∫d4x-gxδx0-T-δx0-t0ϕ¯iacabϕ¯˙ib.In (40), the interactions are included in the first term as well as in the second term. In the action above, we have also taken into account the surface term at the boundary which corresponds to the last term of (41). T and t0 in (41) are the upper bound and the lower bound of the time integration, respectively.


## 3.2. Schwinger Dyson Equations

Now let us derive SDEs for both Green’s function and field. These equations can be obtained by taking the variation of the 2PI EA,Γ2, with respect to the scalar field ϕ¯ and Green’s function G.In the following, we first derive SDEs for the field. The variation of the 2PI EA in (39) with respect to the scalar field ϕ¯ leads to(42)1-gxδΓ2δϕ¯iax=-cabJibx-∫d4zcabKijbcx,zccd-gzϕ¯jdz.Using (25) and (26), one computes the right-hand side of the above equation as(43)cabJibx+∫d4zcabKijbcx,zccd-gzϕ¯jdz=-iδx0-t0vim~itanh⁡βm~i2+cabκiibck=0ccdat03vid=0,where we have used κab(k) given in (23). The left-hand side of Eq. (42) is computed using (40) and one obtains the following equation of motion of the scalar field ϕ¯,(44)δij□+m~ij2ϕ¯jdx=cdaDabcAijkϕ¯jbxϕ¯kcx+Gjkbcx,x,where the Laplacian of Friedman-Lemaître-Robertson-Walker metric is given by(45)□=∇μ∇μ=∂2∂x02-1ax02∇·∇+3a˙a∂∂x0.Next, the equation of motion for Green’s function is derived in the following way. The variation of the 2PI EA in (39) with respect to Green’s function G leads to(46)δΓ2δGijabx,y=-12cac-gxKijcdx,y-gycdb.The left-hand side of the above equation is obtained by taking variation of Eq. (40) with respect to Green’s function as,(47)δΓ2δGijabx,y=-i2G-1jibay,x+12δ2Sϕ¯,gδϕ¯iaxδϕ¯jby,where the second term of above expression is computed using action in (41). Taking all together (46) and (47), one obtains the following two differential equations for Green’s function,(48)□→x+m~i2Gijabx,y=-iδijcab-gxδx-y+2cadDdceAiklϕ¯l,xeGkj,xycb+∫d4zKikaex,z-gzcefGkjfbz,y,(49)□→y+m~j2Gijabx,y=-iδijδx-ycab-gy+2Gik,xyacDcefAkjlϕ¯l,yfceb+∫d4zGikaex,zcef-gzKkjfbz,y,where □x=∇xμ∇μx and □y=∇yμ∇μy.Next, we rescale Green’s function, field and coupling constant of interaction as follows:(50)ϕ¯x0≕at0ax03/2φ^x0,(51)Gx0,y0,k≕at0ax03/2G^x0,y0,kat0ay03/2,(52)A^x0≔at0ax03/2A,where at0 stands for the initial value for the scale factor and we have defined at0≔a(t0) and we have used Fourier transformation for Green’s function as(53)Gx0,y0,k=∫d3rGx0,r,y0,0eik·r.By using these new definitions, SDEs in (44) are written as(54)∂2∂x02+m¯i2x0φ^idx0=cdaDabcA^ijkx0φ^jbx0φ^kcx0+G^jkbcx,x.Next SDEs for the rescaled Green’s function in (48) and (49) are written as(55)∂2∂x02+Ωi,k2x0G^ij,x0y0abk=2cadDdceA^ikl,x0φ^l,x0eG^kj,x0y0cbk-iδijδx0y0cabat03-iδt0x0κikaekat03cefG^kj,t0y0fbk,(56)∂2∂y02+Ωi,k2y0G^ij,x0y0abk=2G^ik,x0y0ackDcefA^kjl,y0φ^l,y0fceb-iδijδx0y0cabat03-iG^ik,x0t0aekcefat03κkjfbkδt0y0,where we have defined(57)Ωi,k2x0≔k2ax02+m¯i2x0,(58)m¯i2x0≔m~i2x0-32a¨x0ax0-34a˙x0ax02.Note that the first derivative with respect to time which is originally presented in the expression of Laplacian (45) is now absent in the expression of SDEs for the rescaled fields and Green’s functions.


## 3.3. The Initial Condition for Green’s Function and Field

In this subsection, the initial conditions for Green’s function and field are determined. For simplicity, let us look back to example model for one real scalar field. We first compute the initial condensation of the fieldϕ¯(t0)=φ^(t0) (see (50)). Using (32) and setting x0=t0, we compute it as follows:(59)ϕ¯t0≡ϕt0,x≔∫dϕϕt0,xexp⁡-SEclexampleϕ,ϕ∫dϕexp⁡-SEclexampleϕ,ϕ=∫dϕϕxexp⁡-1/2∫d3xd3yϕxDx-yϕy+2at03∫d3xϕxj1∫dϕexp⁡-1/2∫d3xd3yϕxDx-yϕy+2at03∫d3xϕxj1,where we have computed the last term of (22) using (25) and D(r) is defined as [14](60)Dr=2at03∫d3k2π3ωkcosh⁡βωk-1sinh⁡βωke-ir·k.To proceed the calculation, we denote J(x) as J(x)=2at03j1. Then the initial condensation of field 〈ϕ(t0,x)〉 is given by(61)ϕt0,x=∫dϕ′xϕ′x+∫d3zD-1x-zJzexp-1/2∫d3xd3yϕ′xDx-yϕ′y∫dϕ′xexp-1/2∫d3xd3yϕ′xDx-yϕ′y=v,where we have defined ϕ′(x)=ϕ(x)-∫d3zD-1(z-x)J(z) and D-1(z-x) satisfies(62)∫d3xD-1z-xDx-y=δz-y.Next we will compute the initial condition for Green’s functionG(t0,t0,k)=G^(t0,t0,k) (see (51)). By using (34) and setting x0=t0 and y0=t0, one computes it as follows:(63)Gabt0,x,t0,y=∫∫dϕ1dϕ2ϕ2Φt0,yΦt0,xϕ1exp⁡-SEclexampleϕ1,ϕ2∫dϕexp⁡-SEclexampleϕ,ϕ=∫dϕ′ϕ′yϕ′xexp⁡-1/2∫d3xd3yϕ′xDx-yϕ′y∫dϕ′exp⁡-1/2∫d3xd3yϕ′xDx-yϕ′y,=D-1y-x.Using (53) and (60), (63) becomes(64)Gabt0,t0,k=D-1k=12ωkat03sinh⁡βωkcosh⁡βωk-1.The above results with the example model can be extended to our model and we summarize them as(65)G^ij,t0t0abk=δij12ωikat03sinh⁡βωikcosh⁡βωik-1,(66)φ^iat0=vi.Next we derive the time derivative of the field and Green’s function at the initial timet0. First we integrate the field equation in (54) with respect to time. By setting x0=t0, we obtain(67)∂φ^i,x0∂x0x0=t0=0.Similarly, we integrate (55) with respect to time x0. By setting both x0 and y0 equal to t0, we obtain the following initial condition(68)limx0→t0⁡∂∂x0G^ij,x0t0abk=-iδijcabat03-iκikaekat03cefG^kj,t0t0fbk.Finally, we integrate (56) with respect to time y0. By setting both x0 and y0 equal to t0, we obtain another initial condition,(69)limy0→t0⁡∂∂y0G^ij,t0y0abk=-iδijcabat03-iG^ik,t0t0aekcefat03κkjfbk.


## 4. The Expectation Value of PNA

The SDEs obtained in the previous section allow us to write the solutions for both Green’s functions and fields in the form of integral equations. In this section, we present the correction to the expectation value of the PNA up to the first order contribution with respect to the cubic interaction. For this purpose, in Subsection4.1, we show how one analytically obtains the solutions of SDEs. We write down the solutions up to the first order of the cubic interaction. In the Subsection 4.2, we also write the expectation value of the PNA up to the first order of the cubic interaction and investigate it by taking into account of the time dependence of the scale factor.


### 4.1. The Solution of Green’s Function and Fields includingo(A) Corrections

The SDEs in present work are inhomogeneous differential equations of the second order. To solve the differential equation, the variation of constants method is used. With the method, the solutions of SDEs are written in the form of integral equations. We solve the integral equation perturbatively and the solutions up to the first order of the cubic interaction are obtained. We first write the solutions of fields as(70)φ^i,x0d=φ^i,x0d,free+φ^i,x0d,oA,(71)φ^i,x0d,free=-K¯i,x0t0′φ^i,t0d,(72)φ^i,x0d,oA=∫t0x0K¯i,x0tcdaDabcA^ijktφ^jb,freetφ^kc,freet+∫d3k2π3G^jk,ttbc,freekdt,where φ^free denotes the free part contribution while φ^o(A) is the contribution due to the first order of the cubic interaction. In Appendix A.2, (70)-(72) are derived in detail. K¯i,x0y0≔K¯i,x0y0,k=0 and K¯i,x0y0,k is defined by(73)K¯i,x0y0,k≔1Wi,kfi,kx0gi,ky0-gi,kx0fi,ky0,where Wi,k is defined as(74)Wi,k≔f˙i,kx0gi,kx0-fi,kx0g˙i,kx0.fi,k and gi,k are the solutions which satisfy the following homogeneous differential equations:(75)∂2∂x02+Ωi,k2x0fi,kx0=0,(76)∂2∂x02+Ωi,k2x0gi,kx0=0,where Ωi,k2(x0) is given in (57). In Appendix B, fi,k and gi,k are derived in detail. K′¯i,x0y0≔K′¯i,x0y0,k=0 and K′¯i,x0y0,k is also defined as follows:(77)K′¯i,x0y0,k≔∂K¯i,x0y0,k∂y0.Next we write down the solution of Green’s function as follows:(78)G^ij,x0y0abk=G^ij,x0y0ab,freek+G^ij,x0y0ab,oAk,(79)G^ij,x0y0ab,freek=δij2ωi,kat03coth⁡βωi,k21111abK¯i,x0t0,k′K¯i,y0t0,k′+ωi,k2K¯i,x0t0,kK¯y0t0,k+iδij2at03K¯i,x0y0,kϵab+cabθy0-x0-θx0-y0,where ϵab is an antisymmetric tensor and its nonzero components are given as ϵ12=1 while θ(t) denotes a unit step function,(80)G^ij,x0y0ab,oAk=∫t0y0Rij,x0t,kab,oAK¯j,y0t,kdt+∫t0x0K¯i,x0t,kQij,tt0,kac,oAEjj,kT,cbK¯j,y0t0,k-K¯j,y0t0,k′δcbdt,where Qo(A), Ro(A), and Ek are given as(81)Qij,x0y0,kab,oA=2cadDdceA^ikl,x0φ^l,x0e,freeG^kj,x0y0cb,freek,(82)Rij,x0y0,kab,oA=2G^ik,x0y0ac,freekDcefA^kjl,y0φ^l,y0f,freeceb,(83)Eik,kac=-iκikaekat03cec,and κijab(k) is given in (23). In Appendices A.4 and A.5, we derive (79) and (80) in detail, respectively.


### 4.2. The Expectation Value of PNA includingo(A) Corrections and the First Order of the Hubble Parameter

Next we compute the PNA in (36) including the first order correction with respect to A (o(A)) and the effect of expansion up to the first order of the Hubble parameter. By using rescaled fields, Green’s function, and coupling constant in (50)-(52), one can write down total contribution to the expectation value of PNA with order o(A) corrections as(84)
a
x
0
a
t
0
3
j
0
x
0
=
R
e
φ
^
2
1
,
f
r
e
e
x
0
φ
^
˙
1
1
∗
,
f
r
e
e
x
0
-
φ
^
1
1
,
f
r
e
e
x
0
φ
^
˙
2
1
∗
,
f
r
e
e
x
0
+
∫
d
3
k
2
π
3
∂
∂
x
0
-
∂
∂
y
0
R
e
G
^
12
12
,
o
A
x
0
,
y
0
,
k
y
0
→
x
0
+
R
e
φ
^
2
1
,
f
r
e
e
x
0
φ
^
˙
1
1
∗
,
o
A
x
0
-
φ
^
1
1
,
f
r
e
e
x
0
φ
^
˙
2
1
∗
,
o
A
x
0
+
R
e
φ
^
2
1
,
o
A
x
0
φ
^
˙
1
1
∗
,
f
r
e
e
x
0
-
φ
^
1
1
,
o
A
x
0
φ
^
˙
2
1
∗
,
f
r
e
e
x
0The first line of the above equation is the zeroth order of the cubic interaction while the next three terms are the first order.As was indicated previously, we will further investigate the expectation value of the PNA for the case of time-dependent scale factor. For that purpose, one can expand scale factor aroundt0 for 0<t0⩽x0 as follows: (85)ax0=at0+x0-t0a˙t0+12x0-t02a¨t0+…=a0+a1x0+a2x0+….We first assume that a(n+1)(x0)<a(n)(x0) when x0 is near t0. Then one can keep only the following terms:(86)ax0≃a0+a1x0,and a(n)(x0) for (n⩾2) are set to be zero. a(0) corresponds to the constant scale factor and a(1)(x0) corresponds to linear Hubble parameter H(t0). Thus it can be written as(87)ax0at0=1+x0-t0Ht0,where H(t0) is given by(88)Ht0=a˙t0at0,and t0>0. Throughout this study, we only keep first order of H(t0) as the first nontrivial approximation. For the case that Hubble parameter is positive, it corresponds to the case for the expanding universe. Under this situation, a˙(x0)=a(t0)H(t0) and a¨(x0)=0.Now let us briefly go back to (58). With these approximations, the second term of (58) is apparently vanished. Since a˙(x0) is proportional to linear H(t0), the third term of (58) involves second order of H(t0). Hence, one can neglect it and the Riemann curvature R(x0) in (3) is also vanished. Therefore, m¯i2(x0) is simply written as m~i2. Now m~i2 are given as(89)m~12=mϕ2-B2,(90)m~22=mϕ2+B2,(91)m~32=mN2.Next we define ωi,k as(92)ωi,k≔k2at02+m~i2.We consider Ωi,k(x0) defined in (57). One can expand it around time t0 as(93)Ωi,kx0≃ωi,k+x0-t0∂∂x0Ωi,kx0x0=t0=ωi,k1-Ht0x0-t0k2at0ωi,kt02.Now let us investigate the expectation value of PNA under these approximations. For the case thatφ^1,t0=φ^2,t0=0 and φ^3,t0≠0, the nonzero contribution to the expectation value of PNA comes only from o(A) corrections to Green’s function. From (84), we can obtain(94)j0x0=2ax03φ^3,t0∫d3k2π3∫t0x0A^123,t-K¯3,tt0,0′12ω2,kt0×coth⁡βω2,kt02K¯˙1,x0t,kK¯2,x0t0,k′-K¯1,x0t,kK¯˙2,x0t0,k′K¯2,tt0,k′+ω2,k2t0K¯˙1,x0t,kK¯2,x0t0,k-K¯1,x0t,kK¯˙2,x0t0,kK¯2,tt0,k-1↔2forlowerindicesdt,where we have used (71) and (80). Following the expression of the scale factor in (87), K¯ is also divided into the part of the constant scale factor and the part which is proportional to H(t0). In Appendix C, K¯ and its derivative are derived in detail. In the above expression, H(t0) is also included in A^(t). Since we are interested in the PNA up to the first order of H(t0), we expand it as follows:(95)A^t≃A1-32t-t0Ht0.Furthermore, substituting (87), (95) and K¯ and its derivative in (C.7), (C.8), (C.11)-(C.16) into Eq. (94), one can divide the PNA into two parts,(96)j0x0=j0x01st+j0x02nd,(97)j0x01st=2φ^3,t0A123at03∫d3k2π3∫t0x01-3x0-t0Ht0-32t-t0Ht0×-K¯3,tt0,00′2ω2,kt0coth⁡βω2,kt02K¯2,x0t0,k0′∂.⃡K¯1,x0t,k0K¯2,tt0,k0′+ω2,k2t0K¯2,x0t0,k0∂.⃡K¯1,x0t,k0K¯2,tt0,k0-1↔2for  lower  indicesdt,(98)j0x02nd=2φ^3,t0A123at03∫d3k2π3∫t0x0-K¯3,tt0,00′2ω2,kt0coth⁡βω2,kt02×K¯2,x0t0,k0′∂.⃡K¯1,x0t,k0K¯2,tt0,k1′+K¯2,x0t0,k1′∂.⃡K¯1,x0t,k0+K¯2,x0t0,k0′∂.⃡K¯1,x0t,k1K¯2,tt0,k0′+ω2,k2t0K¯2,x0t0,k0∂.⃡K¯1,x0t,k0K¯2,tt0,k1+K¯2,x0t0,k1∂.⃡K¯1,x0t,k0+K¯2,x0t0,k0∂.⃡K¯1,x0t,k1K¯2,tt0,k0-1↔2forlowerindicesdt,and the derivative ∂.⃡ acts on the first argument of K¯ and is defined as follows:(99)K¯2,x0t,k∂.⃡K¯1,x0t,k=K¯2,x0t,k∂∂x0K¯1,x0t,k-∂∂x0K¯2,x0t,kK¯1,x0t,k.Each term of the PNA shown in (97) and (98) can be understood as follows. The first term is the PNA with the constant scale factor. The second term with a prefactor -3H(t0)(x0-t0)1/at03≃1/a(x0)3-1/at03 is called the dilution effect. The third term with a prefactor -3/2A123(t-t0)H(t0)≃A^123(t)-A123 is called the freezing interaction effect. The fourth term which corresponds to 〈j0(x0)〉2nd is called the redshift effect. Below we explain their physical origins. The dilution of the PNA is caused by the increase of the volume of the universe. The origin of the freezing interaction effect can be understood with (95). It implies that the strength of the cubic interaction A^(t) controlling the size of PNA decreases as the scale factor grows. The origin of the redshift can be explained as follows. As shown in (57), as the scale factor grows, the physical wavelength becomes large. Therefore, the momentum and the energy of the particles become small. Note that this effect does not apply to the zero-mode such as condensate which is homogeneous and is a constant in the space.Before closing this section, we compute the production rate of the PNA per unit time which is a useful expression when we understand the numerical results of the PNA. We compute the time derivative of the PNA for the case of the constant scale factorHt0=0. By setting Ht0=0, one obtains it at the initial time x1=x0-t0=0,(100)∂∂x1j0x1+t0x1=0=v3A123at03∫0∞k2dk2π21ω1,kω2,k×n2-n1ω1,k+ω2,k+n2+n1+1ω1,k-ω2,k,where ni is the distribution functions for the Bose particles,(101)ni=e-βωi,k1-e-βωi,k,i=1,2.In Appendix D, we derive (100) in detail. Because we assume m~1<m~2, one obtains inequality n2<n1. From the expression above, the production rate of PNA at the initial time is negative for v3A123>0. One also finds the rate is logarithmically divergent for the momentum (k) integration,(102)∂∂x1j0x1+t0x1=0divergent≃v3A1232π2at03m~12-m~222log⁡kmaxμ,where μ=O(m~i) (i=1,2) and kmax is an ultraviolet cut off for the momentum integration. With the expression, one expects that for the positive v3A123, the PNA becomes negative from zero just after the initial time and the behavior will be confirmed in the numerical simulation.


## 4.1. The Solution of Green’s Function and Fields includingo(A) Corrections

The SDEs in present work are inhomogeneous differential equations of the second order. To solve the differential equation, the variation of constants method is used. With the method, the solutions of SDEs are written in the form of integral equations. We solve the integral equation perturbatively and the solutions up to the first order of the cubic interaction are obtained. We first write the solutions of fields as(70)φ^i,x0d=φ^i,x0d,free+φ^i,x0d,oA,(71)φ^i,x0d,free=-K¯i,x0t0′φ^i,t0d,(72)φ^i,x0d,oA=∫t0x0K¯i,x0tcdaDabcA^ijktφ^jb,freetφ^kc,freet+∫d3k2π3G^jk,ttbc,freekdt,where φ^free denotes the free part contribution while φ^o(A) is the contribution due to the first order of the cubic interaction. In Appendix A.2, (70)-(72) are derived in detail. K¯i,x0y0≔K¯i,x0y0,k=0 and K¯i,x0y0,k is defined by(73)K¯i,x0y0,k≔1Wi,kfi,kx0gi,ky0-gi,kx0fi,ky0,where Wi,k is defined as(74)Wi,k≔f˙i,kx0gi,kx0-fi,kx0g˙i,kx0.fi,k and gi,k are the solutions which satisfy the following homogeneous differential equations:(75)∂2∂x02+Ωi,k2x0fi,kx0=0,(76)∂2∂x02+Ωi,k2x0gi,kx0=0,where Ωi,k2(x0) is given in (57). In Appendix B, fi,k and gi,k are derived in detail. K′¯i,x0y0≔K′¯i,x0y0,k=0 and K′¯i,x0y0,k is also defined as follows:(77)K′¯i,x0y0,k≔∂K¯i,x0y0,k∂y0.Next we write down the solution of Green’s function as follows:(78)G^ij,x0y0abk=G^ij,x0y0ab,freek+G^ij,x0y0ab,oAk,(79)G^ij,x0y0ab,freek=δij2ωi,kat03coth⁡βωi,k21111abK¯i,x0t0,k′K¯i,y0t0,k′+ωi,k2K¯i,x0t0,kK¯y0t0,k+iδij2at03K¯i,x0y0,kϵab+cabθy0-x0-θx0-y0,where ϵab is an antisymmetric tensor and its nonzero components are given as ϵ12=1 while θ(t) denotes a unit step function,(80)G^ij,x0y0ab,oAk=∫t0y0Rij,x0t,kab,oAK¯j,y0t,kdt+∫t0x0K¯i,x0t,kQij,tt0,kac,oAEjj,kT,cbK¯j,y0t0,k-K¯j,y0t0,k′δcbdt,where Qo(A), Ro(A), and Ek are given as(81)Qij,x0y0,kab,oA=2cadDdceA^ikl,x0φ^l,x0e,freeG^kj,x0y0cb,freek,(82)Rij,x0y0,kab,oA=2G^ik,x0y0ac,freekDcefA^kjl,y0φ^l,y0f,freeceb,(83)Eik,kac=-iκikaekat03cec,and κijab(k) is given in (23). In Appendices A.4 and A.5, we derive (79) and (80) in detail, respectively.


## 4.2. The Expectation Value of PNA includingo(A) Corrections and the First Order of the Hubble Parameter

Next we compute the PNA in (36) including the first order correction with respect to A (o(A)) and the effect of expansion up to the first order of the Hubble parameter. By using rescaled fields, Green’s function, and coupling constant in (50)-(52), one can write down total contribution to the expectation value of PNA with order o(A) corrections as(84)
a
x
0
a
t
0
3
j
0
x
0
=
R
e
φ
^
2
1
,
f
r
e
e
x
0
φ
^
˙
1
1
∗
,
f
r
e
e
x
0
-
φ
^
1
1
,
f
r
e
e
x
0
φ
^
˙
2
1
∗
,
f
r
e
e
x
0
+
∫
d
3
k
2
π
3
∂
∂
x
0
-
∂
∂
y
0
R
e
G
^
12
12
,
o
A
x
0
,
y
0
,
k
y
0
→
x
0
+
R
e
φ
^
2
1
,
f
r
e
e
x
0
φ
^
˙
1
1
∗
,
o
A
x
0
-
φ
^
1
1
,
f
r
e
e
x
0
φ
^
˙
2
1
∗
,
o
A
x
0
+
R
e
φ
^
2
1
,
o
A
x
0
φ
^
˙
1
1
∗
,
f
r
e
e
x
0
-
φ
^
1
1
,
o
A
x
0
φ
^
˙
2
1
∗
,
f
r
e
e
x
0The first line of the above equation is the zeroth order of the cubic interaction while the next three terms are the first order.As was indicated previously, we will further investigate the expectation value of the PNA for the case of time-dependent scale factor. For that purpose, one can expand scale factor aroundt0 for 0<t0⩽x0 as follows: (85)ax0=at0+x0-t0a˙t0+12x0-t02a¨t0+…=a0+a1x0+a2x0+….We first assume that a(n+1)(x0)<a(n)(x0) when x0 is near t0. Then one can keep only the following terms:(86)ax0≃a0+a1x0,and a(n)(x0) for (n⩾2) are set to be zero. a(0) corresponds to the constant scale factor and a(1)(x0) corresponds to linear Hubble parameter H(t0). Thus it can be written as(87)ax0at0=1+x0-t0Ht0,where H(t0) is given by(88)Ht0=a˙t0at0,and t0>0. Throughout this study, we only keep first order of H(t0) as the first nontrivial approximation. For the case that Hubble parameter is positive, it corresponds to the case for the expanding universe. Under this situation, a˙(x0)=a(t0)H(t0) and a¨(x0)=0.Now let us briefly go back to (58). With these approximations, the second term of (58) is apparently vanished. Since a˙(x0) is proportional to linear H(t0), the third term of (58) involves second order of H(t0). Hence, one can neglect it and the Riemann curvature R(x0) in (3) is also vanished. Therefore, m¯i2(x0) is simply written as m~i2. Now m~i2 are given as(89)m~12=mϕ2-B2,(90)m~22=mϕ2+B2,(91)m~32=mN2.Next we define ωi,k as(92)ωi,k≔k2at02+m~i2.We consider Ωi,k(x0) defined in (57). One can expand it around time t0 as(93)Ωi,kx0≃ωi,k+x0-t0∂∂x0Ωi,kx0x0=t0=ωi,k1-Ht0x0-t0k2at0ωi,kt02.Now let us investigate the expectation value of PNA under these approximations. For the case thatφ^1,t0=φ^2,t0=0 and φ^3,t0≠0, the nonzero contribution to the expectation value of PNA comes only from o(A) corrections to Green’s function. From (84), we can obtain(94)j0x0=2ax03φ^3,t0∫d3k2π3∫t0x0A^123,t-K¯3,tt0,0′12ω2,kt0×coth⁡βω2,kt02K¯˙1,x0t,kK¯2,x0t0,k′-K¯1,x0t,kK¯˙2,x0t0,k′K¯2,tt0,k′+ω2,k2t0K¯˙1,x0t,kK¯2,x0t0,k-K¯1,x0t,kK¯˙2,x0t0,kK¯2,tt0,k-1↔2forlowerindicesdt,where we have used (71) and (80). Following the expression of the scale factor in (87), K¯ is also divided into the part of the constant scale factor and the part which is proportional to H(t0). In Appendix C, K¯ and its derivative are derived in detail. In the above expression, H(t0) is also included in A^(t). Since we are interested in the PNA up to the first order of H(t0), we expand it as follows:(95)A^t≃A1-32t-t0Ht0.Furthermore, substituting (87), (95) and K¯ and its derivative in (C.7), (C.8), (C.11)-(C.16) into Eq. (94), one can divide the PNA into two parts,(96)j0x0=j0x01st+j0x02nd,(97)j0x01st=2φ^3,t0A123at03∫d3k2π3∫t0x01-3x0-t0Ht0-32t-t0Ht0×-K¯3,tt0,00′2ω2,kt0coth⁡βω2,kt02K¯2,x0t0,k0′∂.⃡K¯1,x0t,k0K¯2,tt0,k0′+ω2,k2t0K¯2,x0t0,k0∂.⃡K¯1,x0t,k0K¯2,tt0,k0-1↔2for  lower  indicesdt,(98)j0x02nd=2φ^3,t0A123at03∫d3k2π3∫t0x0-K¯3,tt0,00′2ω2,kt0coth⁡βω2,kt02×K¯2,x0t0,k0′∂.⃡K¯1,x0t,k0K¯2,tt0,k1′+K¯2,x0t0,k1′∂.⃡K¯1,x0t,k0+K¯2,x0t0,k0′∂.⃡K¯1,x0t,k1K¯2,tt0,k0′+ω2,k2t0K¯2,x0t0,k0∂.⃡K¯1,x0t,k0K¯2,tt0,k1+K¯2,x0t0,k1∂.⃡K¯1,x0t,k0+K¯2,x0t0,k0∂.⃡K¯1,x0t,k1K¯2,tt0,k0-1↔2forlowerindicesdt,and the derivative ∂.⃡ acts on the first argument of K¯ and is defined as follows:(99)K¯2,x0t,k∂.⃡K¯1,x0t,k=K¯2,x0t,k∂∂x0K¯1,x0t,k-∂∂x0K¯2,x0t,kK¯1,x0t,k.Each term of the PNA shown in (97) and (98) can be understood as follows. The first term is the PNA with the constant scale factor. The second term with a prefactor -3H(t0)(x0-t0)1/at03≃1/a(x0)3-1/at03 is called the dilution effect. The third term with a prefactor -3/2A123(t-t0)H(t0)≃A^123(t)-A123 is called the freezing interaction effect. The fourth term which corresponds to 〈j0(x0)〉2nd is called the redshift effect. Below we explain their physical origins. The dilution of the PNA is caused by the increase of the volume of the universe. The origin of the freezing interaction effect can be understood with (95). It implies that the strength of the cubic interaction A^(t) controlling the size of PNA decreases as the scale factor grows. The origin of the redshift can be explained as follows. As shown in (57), as the scale factor grows, the physical wavelength becomes large. Therefore, the momentum and the energy of the particles become small. Note that this effect does not apply to the zero-mode such as condensate which is homogeneous and is a constant in the space.Before closing this section, we compute the production rate of the PNA per unit time which is a useful expression when we understand the numerical results of the PNA. We compute the time derivative of the PNA for the case of the constant scale factorHt0=0. By setting Ht0=0, one obtains it at the initial time x1=x0-t0=0,(100)∂∂x1j0x1+t0x1=0=v3A123at03∫0∞k2dk2π21ω1,kω2,k×n2-n1ω1,k+ω2,k+n2+n1+1ω1,k-ω2,k,where ni is the distribution functions for the Bose particles,(101)ni=e-βωi,k1-e-βωi,k,i=1,2.In Appendix D, we derive (100) in detail. Because we assume m~1<m~2, one obtains inequality n2<n1. From the expression above, the production rate of PNA at the initial time is negative for v3A123>0. One also finds the rate is logarithmically divergent for the momentum (k) integration,(102)∂∂x1j0x1+t0x1=0divergent≃v3A1232π2at03m~12-m~222log⁡kmaxμ,where μ=O(m~i) (i=1,2) and kmax is an ultraviolet cut off for the momentum integration. With the expression, one expects that for the positive v3A123, the PNA becomes negative from zero just after the initial time and the behavior will be confirmed in the numerical simulation.


## 5. Numerical Results

In this section, we numerically study the time dependence of the PNA. The PNA depends on the parameters of the model such as masses and coupling constants. It also depends on the initial conditions and the expansion rates of the universe. Since the PNA is linearly proportional to the coupling constantA123 and the initial value of the field φ^3,t0, we can set these parameters as unity in the unit of energy and later on one can multiply their values. As for the initial scale factor at0, without loss of generality, one can set this dimensionless factor as unity. For the other parameters of the model, we choose m~2,B, and ω3,0=m~3 as independent parameters since the mass m~1 is written as(103)m~12=m~22-2B2.The temperature T and the expansion rate H(t0) determine the environment for the universe. The former determines the thermal distribution of the scalar fields. Within the approximation for the time dependence of the scale factor in (87), H(t0) is the only parameter which controls the expansion rate of the universe. The approximation is good for the time range which satisfies the following inequality:(104)x0-t0≪13Ht0.The time dependence of PNA is plotted as a function of the dimensionless time defined as(105)t=ω3,0rx0-t0,where ω3,0r is a reference frequency. In terms of the dimensionless time, the condition of (104) is written as(106)t≪tmax≡ω3,0r3Ht0.How the PNA behaves with respect to time is discussed in the following Subsections5.1–5.3. The results, as will be shown later, revealed that the PNA has an oscillatory behavior. We also investigate the parameter dependence for two typical cases, one of which corresponds to the longer period and the other corresponds to the shorter period. In the numerical simulation, we do not specify the unit of parameters. Note that the numerical values for the dimensionless quantities such as ratio of masses do not depend on the choice of the unit as far as the quantities in the ratio are given in the same unit. In Subsection 5.4, we assign the unit for the parameters and estimate the ratio of the PNA over entropy density.


### 5.1. The PNA with the Longer Period

Let us now consider the PNA which has the longer period. While we investigate the dependence of several parameters, we fix two parameters as(m~2,Ht0)=(0.05,10-3). In Figure 1, the temperature (T) dependence of PNA is shown. It depends on the temperature only through hyperbolic function as shown in (94). In this figure, tmax in (106) is around 110. What stands out of this figure is the change of the amplitude for PNA among the three curves. As the temperature increases, the amplitude of the oscillation becomes larger. In Figure 2, we show the B dependence. Interestingly, both of the amplitude and the period of the oscillation change when we alter the parameter B. As it increases, the amplitude becomes larger and its period becomes shorter. Figure 3 shows the dependence of the PNA on ω3,0. As shown in the black, blue, and dot-dashed blue lines, the position of the first node does not change when ω3,0 takes its value within the difference of m~1 and m~2. However, the amplitude of oscillation gradually decreases as ω3,0 increases up to the mass difference. The more interesting findings were observed when ω3,0 becomes larger than the mass difference. As ω3,0 becomes larger, the amplitude decreases and the new node is formed at once. The dashed and dotted blue lines show this behavior. The dependence on the expansion rate (Ht0) is shown in Figure 4. There is an interesting aspect of this figure at the fixed time t. As the expansion rate becomes larger, the size of PNA becomes smaller.Figure 1
Dependence on temperatureT of the time evolution of PNA. In horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,Ht0,ω3,0) = (0.04,0.05,0.021,10-3,0.0035) for all of the lines. The dotted red, black, dashed red, and red lines show the cases T=50,100,200, and 400, respectively.Figure 2
Dependence on parameterB of the time evolution of PNA. The horizontal axis is the dimensionless time defined as t=ω3,0r(x0-t0). As a reference angular frequency, we choose ω3,0r=0.35. We fix a set of parameters as (m~2,T,Ht0,ω3,0) = (0.05,100,10-3,0.0035) for all of the lines. The purple, thin purple, dotted purple, black, dashed purple, and dot-dashed purple lines show the cases B=0.027,0.025,0.024,0.021,0.018 and 0.012, respectively.Figure 3
Theω3,0 dependence of the time evolution of PNA. In horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We use a set of parameters as (m~1,m~2,B,T,Ht0) = (0.04,0.05,0.021,100,10-3) for all of the lines. The black, blue, dot-dashed blue, dashed blue and dotted blue lines show the cases ω3,0=0.0035,0.0045,0.008,0.012 and 0.015, respectively.Figure 4
Dependence on the expansion rateHt0 of the time evolution of PNA. In the horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,T,ω3,0) = (0.04,0.05,0.021,100,0.0035) for all of the lines. The gray, dashed gray, dotted gray, and black lines show the cases Ht0=0,10-4,5×10-4 and 10-3, respectively.


### 5.2. The PNA with the Shorter Period

Now we investigate the PNA with the shorter period. In Figure5, we show the temperature (T) dependence for the time evolution of PNA. In this regard, the temperature dependence is similar to the one with the longer period. Namely, the amplitude of oscillation becomes larger as the temperature increases. Figure 6 shows the B dependence. As B parameter decreases, the period of oscillation becomes longer. However, there were different effects on the amplitude of oscillation. In the left plot, we show the cases that the mass difference m~2-m~1 is larger than the frequency ω3,0. Since B2, proportional to mass squared difference m~22-m~12, of the magenta line is smaller than that of the black line, the mass difference m~2-m~1 of the magenta line is closer to ω3,0. At the beginning (0<t<22), the black line of large B has the larger amplitude than that of the magenta line of small B. At time t~22, the amplitude of the magenta line becomes larger than that of the black line. We also observed that when the mass difference is near to the ω3,0, that is for the case of magenta line, the amplitude grows slowly compared with that of the black line and reaches its maximal value between one and a half period and twice of the period. After taking its maximal value, it slowly decreases. In the right plot, the blue line shows the case that the mass difference m~2-m~1 is smaller than the frequency ω3,0. In comparison with the black line, the phase shift of π/2 was observed in the blue line. The dependence on the parameter B is similar to that of the magenta line. Namely, as B becomes smaller, the amplitude gradually grows at the beginning and slowly decreases at the later time.Figure 5
Dependence on temperatureT of the time evolution of PNA. As for horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,ω3,Ht0)=(2,3,1.58,0.35,10-3) for all the lines. The black, light red, and red lines show the cases T=100,200, and 400, respectively.Figure 6
B dependence for the time evolution of PNA. The horizontal axis is the dimensionless time defined as t=ω3,0r(x0-t0). As a reference angular frequency, we choose ω3,0r=0.35. We use a set of parameters as (m~2,T,Ht0,ω3,0)=(3,100,10-3,0.35) for all the lines. In the left plot, the black and magenta lines display the cases B=1.58 and 1.01, respectively. For the right plot, the black and blue lines display the cases B=1.58 and 0.76, respectively.In Figure7, we show the dependence on ω3,0. In the left plot, we show the cases that ω3,0’s are smaller than the mass difference as, ω3,0black<ω3,0orange<m~2-m~1. As ω3,0 increases, the period of the oscillation becomes shorter. There is also a different behavior of the amplitudes as follows. At the beginning, the amplitudes of both black and orange lines increase. After that, in comparison with the black lines, the amplitude of the orange line slowly decreases. In the right plot, the green line shows the case that ω3,0 is larger than the mass difference. We observe that the amplitude of the green line is smaller than that of the black line and the period of the green one is shorter than that of the black one. Figure 8 shows the dependence of expansion rate (Ht0). In this plot, the PNA gradually decreases as the expansion rate increases.Figure 7
Dependence on frequencyω3,0 of the time evolution of the PNA. We use the dimensionless time t=ω3,0r(x0-t0) for horizontal axis and as a reference angular frequency, we choose ω3,0r=0.35. We fix parameters (m~1,m~2,B,T,Ht0)=(2,3,1.58,100,10-3) for all the lines. The black and orange lines show the cases ω3,0=0.35 and 0.9, respectively (left plot). The black and green lines show the cases ω3,0=0.35 and 1.5, respectively (right plot).Figure 8
The expansion rateH(t0) dependence of the time evolution PNA. In horizontal axes, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix parameters as (m~1,m~2,B,ω3,0,T)=(2,3,1.58,0.35,100) for all the lines. The light gray, brown, and black lines display the cases H(t0)=0,5×10-4 and 10-3, respectively.


### 5.3. The Comparison of Two Different Periods

In this subsection, we present a comparison of two different periods of the time evolution of the PNA. In Figure9, the black line shows the case of the shorter period and the dotted black line shows the case of the longer one. As can be seen in this figure, the PNA with the shorter period frequently changes the sign and the magnitude also strongly depends on the time. In contrast to the shorter period case, both the sign and magnitude of the longer period case are stable if restricted to the time range t=30~60 in Figure 9.Figure 9
Comparison between two different periods of the time evolution of the PNA. We use the dimensionless timet=ω3,0r(x0-t0) for horizontal axis and as a reference we choose ω3,0r=0.35. We fix parameters (T,Ht0)=(100,10-3) for all the lines. The black (shorter period) and black dotted (longer period) lines show the set parameters (m~1,m~2,B,ω3,0) as (2,3,1.58,0.35) and (0.04,0.05,0.021,0.0035), respectively. Notice that our approximation will break down after t=80.


### 5.4. The Evolution of the PNA with the Scale Factor of a Specific Time Dependence

In this subsection, we interpret the numerical simulation in a specific situation. We assume that the time dependence of the scale factor is given by the one in radiation dominated era. We also specify the unit of the parameters, time and temperature. By doing so, we can clarify implication of the numerical simulation in a more concrete situation.Specifically, the time dependence of the scale factor is given as follows:(107)ax0=1+2Ht0x0-t0.The above equation is derived as follows. The Einstein’s equations without cosmological constant lead to the following equation:(108)a˙a2=8π3Gρ,where G is the Newton’s constant. ρ is the energy density for radiation and it is given by(109)ρx0=ρ0a-4x0,where ρ0 is the initial energy density and we set at0=1. By setting x0=t0 in (108), the initial Hubble parameter is given by(110)Ht02=8π3Gρ0.Then using (110), (108) becomes(111)ddx0ax02=2Ht0.Solving the equation above, one can obtain (107).From the expression in (107), one needs to specify the unit of the Hubble parameter at t0. Through (110), it is related to the initial energy density ρ0. Assuming ρ0 is given by radiation with an effective degree of freedom g∗ and a temperature T(t0), one can write ρ0 as follows:(112)ρ0=g∗π230T4t0.Hereafter, we assume that the temperature of the radiation T(t0) is equal to the temperature T in the density operator for the scalar fields. Then one can write the ratio of the initial Hubble parameter and temperature T as follows:(113)Ht0T=π34πg∗5TMPl,where MPl is the Planck mass and MPl=1.2×1019 (GeV). Then one can write the temperature T in GeV unit as follows:(114)TGeV=3π54πg∗Ht0TMPlGeV.In the numerical simulation, the ratio Ht0/T is given. Therefore, for the given ratio and g∗, the temperature T in terms of GeV unit is determined. Then Ht0 in GeV unit also becomes(115)Ht0GeV=TGeV×Ht0T.The masses of the scalar fields m~i (i=1,2,3) can be also expressed in GeV unit as(116)m~iGeV=Ht0GeV×m~iHt0,where we use the ratios m~i/Ht0 given in the numerical simulation.As an example, we study the implication of the numerical simulation shown in Figure9 by specifying the mass parameter in GeV unit. We also determine the unit of time scale. We first determine the temperature in GeV unit using (114). As for the degree of freedom, we can take g∗≃100 which corresponds to the case that all the standard model particles are regarded as radiation. Then, substituting the ratio Ht0/T=10-5 adapted in Figure 9 to (114) and (115), one obtains T~1013 (GeV) and H(t0)~108 (GeV), respectively. The mass parameters are different between the longer period case (the dotted line) and the shorter period case (the solid line). They are also given in GeV unit shown in Table 2. The time scale Δt=100 corresponds to 3×10-9(GeV)-1 which is about 2×10-33 (sec).Table 2
The mass parameters in GeV unit for both longer and shorter period cases in Figure9.
Mass parameter (GeV) The shorter period The longer period m1 2 × 1 0 11 4 × 1 0 9  m2 3 × 1 0 11 5 × 1 0 9  ω3,0 3.5 × 1 0 10 3.5 × 1 0 8One can also estimate the size of PNA. Here, we consider the maximum value of the PNA for the longer period case in Figure9. We evaluate the ratio of the PNA over entropy density s,(117)j0t≃50s=-5×1011GeVTGeV×A123GeVTGeV×v3GeVTGeV×452π2g∗,where s is given by(118)s=g∗2π245T3GeV3.In (117), the first numerical factor, -5×1011(GeV), is obtained in the following way:(119)A123GeVTGeV=A123T,(120)v3GeVTGeV=v3T,(121)-5×1011GeVTGeV=-5T.In the right-hand side of the above equations, T denotes the temperature in the universal unit of the simulation in Figure 9 where T=100 is used. In the left-hand side, T denotes the corresponding temperature in GeV unit and it is T=1013 (GeV). Substituting the temperature T=1013 (GeV) into (117), one obtains(122)j0t≃50s=-1×10-11×A123GeV108GeV×v3GeV1010GeV.From the equation above, we can achieve the ratio as 10-10 by taking A123=108 (GeV) and v3=1011 (GeV).


## 5.1. The PNA with the Longer Period

Let us now consider the PNA which has the longer period. While we investigate the dependence of several parameters, we fix two parameters as(m~2,Ht0)=(0.05,10-3). In Figure 1, the temperature (T) dependence of PNA is shown. It depends on the temperature only through hyperbolic function as shown in (94). In this figure, tmax in (106) is around 110. What stands out of this figure is the change of the amplitude for PNA among the three curves. As the temperature increases, the amplitude of the oscillation becomes larger. In Figure 2, we show the B dependence. Interestingly, both of the amplitude and the period of the oscillation change when we alter the parameter B. As it increases, the amplitude becomes larger and its period becomes shorter. Figure 3 shows the dependence of the PNA on ω3,0. As shown in the black, blue, and dot-dashed blue lines, the position of the first node does not change when ω3,0 takes its value within the difference of m~1 and m~2. However, the amplitude of oscillation gradually decreases as ω3,0 increases up to the mass difference. The more interesting findings were observed when ω3,0 becomes larger than the mass difference. As ω3,0 becomes larger, the amplitude decreases and the new node is formed at once. The dashed and dotted blue lines show this behavior. The dependence on the expansion rate (Ht0) is shown in Figure 4. There is an interesting aspect of this figure at the fixed time t. As the expansion rate becomes larger, the size of PNA becomes smaller.Figure 1
Dependence on temperatureT of the time evolution of PNA. In horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,Ht0,ω3,0) = (0.04,0.05,0.021,10-3,0.0035) for all of the lines. The dotted red, black, dashed red, and red lines show the cases T=50,100,200, and 400, respectively.Figure 2
Dependence on parameterB of the time evolution of PNA. The horizontal axis is the dimensionless time defined as t=ω3,0r(x0-t0). As a reference angular frequency, we choose ω3,0r=0.35. We fix a set of parameters as (m~2,T,Ht0,ω3,0) = (0.05,100,10-3,0.0035) for all of the lines. The purple, thin purple, dotted purple, black, dashed purple, and dot-dashed purple lines show the cases B=0.027,0.025,0.024,0.021,0.018 and 0.012, respectively.Figure 3
Theω3,0 dependence of the time evolution of PNA. In horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We use a set of parameters as (m~1,m~2,B,T,Ht0) = (0.04,0.05,0.021,100,10-3) for all of the lines. The black, blue, dot-dashed blue, dashed blue and dotted blue lines show the cases ω3,0=0.0035,0.0045,0.008,0.012 and 0.015, respectively.Figure 4
Dependence on the expansion rateHt0 of the time evolution of PNA. In the horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,T,ω3,0) = (0.04,0.05,0.021,100,0.0035) for all of the lines. The gray, dashed gray, dotted gray, and black lines show the cases Ht0=0,10-4,5×10-4 and 10-3, respectively.


## 5.2. The PNA with the Shorter Period

Now we investigate the PNA with the shorter period. In Figure5, we show the temperature (T) dependence for the time evolution of PNA. In this regard, the temperature dependence is similar to the one with the longer period. Namely, the amplitude of oscillation becomes larger as the temperature increases. Figure 6 shows the B dependence. As B parameter decreases, the period of oscillation becomes longer. However, there were different effects on the amplitude of oscillation. In the left plot, we show the cases that the mass difference m~2-m~1 is larger than the frequency ω3,0. Since B2, proportional to mass squared difference m~22-m~12, of the magenta line is smaller than that of the black line, the mass difference m~2-m~1 of the magenta line is closer to ω3,0. At the beginning (0<t<22), the black line of large B has the larger amplitude than that of the magenta line of small B. At time t~22, the amplitude of the magenta line becomes larger than that of the black line. We also observed that when the mass difference is near to the ω3,0, that is for the case of magenta line, the amplitude grows slowly compared with that of the black line and reaches its maximal value between one and a half period and twice of the period. After taking its maximal value, it slowly decreases. In the right plot, the blue line shows the case that the mass difference m~2-m~1 is smaller than the frequency ω3,0. In comparison with the black line, the phase shift of π/2 was observed in the blue line. The dependence on the parameter B is similar to that of the magenta line. Namely, as B becomes smaller, the amplitude gradually grows at the beginning and slowly decreases at the later time.Figure 5
Dependence on temperatureT of the time evolution of PNA. As for horizontal axis, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix a set of parameters as (m~1,m~2,B,ω3,Ht0)=(2,3,1.58,0.35,10-3) for all the lines. The black, light red, and red lines show the cases T=100,200, and 400, respectively.Figure 6
B dependence for the time evolution of PNA. The horizontal axis is the dimensionless time defined as t=ω3,0r(x0-t0). As a reference angular frequency, we choose ω3,0r=0.35. We use a set of parameters as (m~2,T,Ht0,ω3,0)=(3,100,10-3,0.35) for all the lines. In the left plot, the black and magenta lines display the cases B=1.58 and 1.01, respectively. For the right plot, the black and blue lines display the cases B=1.58 and 0.76, respectively.In Figure7, we show the dependence on ω3,0. In the left plot, we show the cases that ω3,0’s are smaller than the mass difference as, ω3,0black<ω3,0orange<m~2-m~1. As ω3,0 increases, the period of the oscillation becomes shorter. There is also a different behavior of the amplitudes as follows. At the beginning, the amplitudes of both black and orange lines increase. After that, in comparison with the black lines, the amplitude of the orange line slowly decreases. In the right plot, the green line shows the case that ω3,0 is larger than the mass difference. We observe that the amplitude of the green line is smaller than that of the black line and the period of the green one is shorter than that of the black one. Figure 8 shows the dependence of expansion rate (Ht0). In this plot, the PNA gradually decreases as the expansion rate increases.Figure 7
Dependence on frequencyω3,0 of the time evolution of the PNA. We use the dimensionless time t=ω3,0r(x0-t0) for horizontal axis and as a reference angular frequency, we choose ω3,0r=0.35. We fix parameters (m~1,m~2,B,T,Ht0)=(2,3,1.58,100,10-3) for all the lines. The black and orange lines show the cases ω3,0=0.35 and 0.9, respectively (left plot). The black and green lines show the cases ω3,0=0.35 and 1.5, respectively (right plot).Figure 8
The expansion rateH(t0) dependence of the time evolution PNA. In horizontal axes, we use the dimensionless time t=ω3,0r(x0-t0) where we choose ω3,0r=0.35. We fix parameters as (m~1,m~2,B,ω3,0,T)=(2,3,1.58,0.35,100) for all the lines. The light gray, brown, and black lines display the cases H(t0)=0,5×10-4 and 10-3, respectively.


## 5.3. The Comparison of Two Different Periods

In this subsection, we present a comparison of two different periods of the time evolution of the PNA. In Figure9, the black line shows the case of the shorter period and the dotted black line shows the case of the longer one. As can be seen in this figure, the PNA with the shorter period frequently changes the sign and the magnitude also strongly depends on the time. In contrast to the shorter period case, both the sign and magnitude of the longer period case are stable if restricted to the time range t=30~60 in Figure 9.Figure 9
Comparison between two different periods of the time evolution of the PNA. We use the dimensionless timet=ω3,0r(x0-t0) for horizontal axis and as a reference we choose ω3,0r=0.35. We fix parameters (T,Ht0)=(100,10-3) for all the lines. The black (shorter period) and black dotted (longer period) lines show the set parameters (m~1,m~2,B,ω3,0) as (2,3,1.58,0.35) and (0.04,0.05,0.021,0.0035), respectively. Notice that our approximation will break down after t=80.


## 5.4. The Evolution of the PNA with the Scale Factor of a Specific Time Dependence

In this subsection, we interpret the numerical simulation in a specific situation. We assume that the time dependence of the scale factor is given by the one in radiation dominated era. We also specify the unit of the parameters, time and temperature. By doing so, we can clarify implication of the numerical simulation in a more concrete situation.Specifically, the time dependence of the scale factor is given as follows:(107)ax0=1+2Ht0x0-t0.The above equation is derived as follows. The Einstein’s equations without cosmological constant lead to the following equation:(108)a˙a2=8π3Gρ,where G is the Newton’s constant. ρ is the energy density for radiation and it is given by(109)ρx0=ρ0a-4x0,where ρ0 is the initial energy density and we set at0=1. By setting x0=t0 in (108), the initial Hubble parameter is given by(110)Ht02=8π3Gρ0.Then using (110), (108) becomes(111)ddx0ax02=2Ht0.Solving the equation above, one can obtain (107).From the expression in (107), one needs to specify the unit of the Hubble parameter at t0. Through (110), it is related to the initial energy density ρ0. Assuming ρ0 is given by radiation with an effective degree of freedom g∗ and a temperature T(t0), one can write ρ0 as follows:(112)ρ0=g∗π230T4t0.Hereafter, we assume that the temperature of the radiation T(t0) is equal to the temperature T in the density operator for the scalar fields. Then one can write the ratio of the initial Hubble parameter and temperature T as follows:(113)Ht0T=π34πg∗5TMPl,where MPl is the Planck mass and MPl=1.2×1019 (GeV). Then one can write the temperature T in GeV unit as follows:(114)TGeV=3π54πg∗Ht0TMPlGeV.In the numerical simulation, the ratio Ht0/T is given. Therefore, for the given ratio and g∗, the temperature T in terms of GeV unit is determined. Then Ht0 in GeV unit also becomes(115)Ht0GeV=TGeV×Ht0T.The masses of the scalar fields m~i (i=1,2,3) can be also expressed in GeV unit as(116)m~iGeV=Ht0GeV×m~iHt0,where we use the ratios m~i/Ht0 given in the numerical simulation.As an example, we study the implication of the numerical simulation shown in Figure9 by specifying the mass parameter in GeV unit. We also determine the unit of time scale. We first determine the temperature in GeV unit using (114). As for the degree of freedom, we can take g∗≃100 which corresponds to the case that all the standard model particles are regarded as radiation. Then, substituting the ratio Ht0/T=10-5 adapted in Figure 9 to (114) and (115), one obtains T~1013 (GeV) and H(t0)~108 (GeV), respectively. The mass parameters are different between the longer period case (the dotted line) and the shorter period case (the solid line). They are also given in GeV unit shown in Table 2. The time scale Δt=100 corresponds to 3×10-9(GeV)-1 which is about 2×10-33 (sec).Table 2
The mass parameters in GeV unit for both longer and shorter period cases in Figure9.
Mass parameter (GeV) The shorter period The longer period m1 2 × 1 0 11 4 × 1 0 9  m2 3 × 1 0 11 5 × 1 0 9  ω3,0 3.5 × 1 0 10 3.5 × 1 0 8One can also estimate the size of PNA. Here, we consider the maximum value of the PNA for the longer period case in Figure9. We evaluate the ratio of the PNA over entropy density s,(117)j0t≃50s=-5×1011GeVTGeV×A123GeVTGeV×v3GeVTGeV×452π2g∗,where s is given by(118)s=g∗2π245T3GeV3.In (117), the first numerical factor, -5×1011(GeV), is obtained in the following way:(119)A123GeVTGeV=A123T,(120)v3GeVTGeV=v3T,(121)-5×1011GeVTGeV=-5T.In the right-hand side of the above equations, T denotes the temperature in the universal unit of the simulation in Figure 9 where T=100 is used. In the left-hand side, T denotes the corresponding temperature in GeV unit and it is T=1013 (GeV). Substituting the temperature T=1013 (GeV) into (117), one obtains(122)j0t≃50s=-1×10-11×A123GeV108GeV×v3GeV1010GeV.From the equation above, we can achieve the ratio as 10-10 by taking A123=108 (GeV) and v3=1011 (GeV).


## 6. Discussion and Conclusion

In this paper, we developed a new mechanism for generating the PNA. This mechanism is realized with the specific model Lagrangian which we have proposed. The model includes a complex scalar. The PNA is associated withU(1) charge of the complex scalar. In addition, we introduce a neutral scalar which interacts with the complex scalar. The U(1) charge is not conserved due to particle number violating interaction. As another source of particle number violation, the U(1) symmetry breaking mass term for the complex scalar is introduced. The initial value for the condensation of the neutral scalar is nonzero. Using 2PI formalism and specifying the initial condition with density operator, the time-dependent PNA is obtained. To include the effect of the time dependence of the scale factor, we approximate it up to the first order of Hubble parameter.The results show that the PNA depends on the interaction couplingA123 and the initial value of the condensation of the neutral scalar φ^3,t0. It also depends on the mass squared difference of two real scalars which originally form a complex scalar. We found that the interaction coupling A123 and the mass squared difference play a key role to give rise to nonvanishing PNA. Even if the initial value of the neutral scalar is nonzero, in the vanishing limit of interaction terms and the mass squared difference, the PNA will vanish. Another important finding is that the contribution to the PNA is divided into four types. The constant scale factor which is the zeroth order of Hubble parameter is the leading contribution. The rest which are the first order term contribute according to their origins. Those are summarized in Table 3.Table 3
The classification ofo(Ht0) contributions to the PNA.
The effect The origin Dilution The increase of volume of the universe due to expansion,1/ax03-1/at03  Freezing interaction The decrease of the strength of the cubic interactionA^ as A^123-A123  Redshift The effective energy of particle as indicated in Eq. (57), k2/a(x0)2+m-i2(x0)We have numerically calculated time evolution of the PNA and have investigated its dependence on the temperature, parameterB, the angular frequency ω3,0, and the expansion rate of the universe. Starting with the null PNA at the initial time, it is generated by particle number violating interaction. Once the nonzero PNA is generated, it starts to oscillate. The amplitude decreases as the time gets larger. The dumping rate of the amplitude increases as the Hubble parameter becomes larger. The period of the oscillation depends on the angular frequency ω3,0 and the parameter B. The former determines the oscillation period for the condensation of the neutral scalar. The latter determines the mass difference of ϕ1 and ϕ2. In the simulation, we focus on the two cases for the oscillation period, one of which corresponds to the longer period case and the other is the shorter period case. The longer period is about half of the Hubble time (1/Ht0) and the shorter period is one percent of the Hubble time. The set of parameters (ω3,0,B) which corresponds to the longer period is typically one percent of the values for the shorter period. In both cases, the amplitude gets larger as the temperature increases. For the longer period, as parameter B becomes larger, the amplitude increases. For the shorter period, in order to have large amplitude, the parameter B is taken so that the mass difference m~2-m~1 is near to ω3,0. In other words, when the resonance condition ω3,0≃m~2-m~1 is satisfied, the amplitude becomes large. For the longer period case, as the angular frequency ω3,0 decreases, the amplitude becomes large.To show how the mechanism can be applied to a realistic situation, we study the simulated results for radiation dominated era when the degree of freedom of light particles is assumed to beg∗≃O(100). Then when the initial temperature of the scalar fields is the same as that of the light particles, the simulation with Ht0/T=10-5 corresponds to the case that the temperature of the universe is 1013 (GeV) which is slightly lower than GUT scale ~1016 (GeV) [20, 21]. The masses of the scalar fields in Figure 9 are different between the shorter period case and the longer period case as shown in Table 2. In the shorter period case, the mass spectrum of the scalar ranges from 1010 (GeV) to 1011 (GeV) while for the longer period case, it is lower than that of the shorter period case by two orders of magnitude. For the longer period case, the maximum asymmetry is achieved at 10-33 (sec) after the initial time. For shorter period, it is achieved at about 10-34 (sec). We have estimated the ratio of the PNA over entropy density by substituting the numerical values of the coupling constant (A123) and the initial expectation value (v3).Compared with the previous works [13, 14, 22, 23], instead of assuming the nonzero PNA at the initial time, the PNA is created through interactions. These interactions have the following unique feature; namely, the interaction between the complex scalars and oscillating condensation of a neutral scalar leads to the PNA. In our work, by assuming the initial condensation of the neutral scalar is away from the equilibrium point, the condensation starts to oscillate. In the expression of the amplitude of PNA, one finds that it is proportional to the CP violating coupling between the scalars and the condensation, the initial condensation of the neutral scalars, and mass difference between mass eigenstates of the two neutral scalars which are originally introduced as a complex scalar with the particle number violating mass and curvature terms. One of the distinctive features of the present mechanism from the one which utilizes the PNA created through the heavy particle decays is as follows. In the mechanism which utilizes the heavy particle decays, the temperature must be high enough so that it once brings the heavy particle to the state of the thermal equilibrium. Therefore the temperature of the universe at reheating era must be as high as the mass of the heavy particle. In contrast to this class of the models, the present model is not restricted by such condition. In place of the condition, the initial condensation must be large enough to explain the asymmetry.In our model, even for the longer period case, the oscillation period is shorter than the Hubble time(1/Ht0). It implies that one of Sakharov conditions for BAU, namely, nonequilibrium condition, is not satisfied. In this respect, we expect that due to the finite life time of the condensation of the neutral scalar, the interaction between the condensation and complex scalars will vanish and eventually the oscillation of the PNA may terminate. The detailed study will be given in the future work. The relation between the PNA and the observed BAU should be also studied. In particular, we need to consider the mechanism how the created PNA is transferred to the observed BAU.



---

*Source: 6825104-2019-04-07.xml*