A new mechanism for generating particle number asymmetry through interactions

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 carries U(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 non-zero. 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.


I. INTRODUCTION
The origin of BAU has long been a question of great interest in explaining why there is more baryon than anti-baryon in nature. Big bang nucleosynthesis (BBN) [1] and cosmic microwave background [2] measurements give the BAU as η ≡ n B /s ⋍ 10 −10 , where n B 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][4][5][6][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][9][10]. Recently, the variety of the method for the calculation of BAU has been also developed [11][12][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 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 non-zero. 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. Especially, the depen-dence 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 section II, we introduce our model with CP and particle number violating interactions. We also specify the density operator as the initial state. In section III, 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 IV, using the solution of Green's function and field, we compute the expectation value of the PNA. Section V provides the numerical study of the time dependence of the PNA. We will also discuss the dependence on the parameters of the model. Section VI 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 non-vanishing Hubble parameter case are derived.

II. A MODEL WITH CP AND PARTICLE NUMBER VIOLATING INTERAC-TION
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 consists of a neutral scalar, N, and a complex scalar, φ. The action is given by, 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 B 2 . 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 N n (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, where a(x 0 ) is the scale factor at time x 0 . Correspondingly the Riemann curvature is given by, In Eq.(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, With these definitions, the free part of the Lagrangian is rewritten as, where the kinetic term is given by,

Cubic interaction coupling Property
and their effective masses are given as follows, Non-zero B 2 or α 2 leads to the non-degenerate mass spectrum for φ 1 and φ 2 . The interaction Lagrangian is rewritten with a totally symmetric coefficient A ijk , with i, j, k = 1, 2, 3. The non-zero components of A ijk are written with the couplings for cubic interaction, A and A 0 , as shown in Table I. We also summarize the qubic interactions and their properties according to U(1) symmetry and CP symmetry.
Nöether current related to the U(1) transformation is, In terms of real scalar fields, the Nöether current alters into, The ordering of the operators in Eq.(11) is arranged so that it is Hermite and the particle number operator, 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, j 0 (x) is the PNA per unit comoving volume. The expectation value of the PNA is written with a density operator, Note that, the PNA is a Heisenberg operator and ρ(t 0 ) is a density operator which specifies the state at the initial time x 0 = t 0 . In this work, we use the density operator with zero chemical potential. It is specifically given by, where β denotes inverse temperature, 1/T , and H is a Hamiltonian which includes linear term of fields, where v i is a constant. The linear term of fields in Eq.(16) is prepared for the non-zero expectation value of fields. Note that the density operator in Eq. (15) is not exactly the same as the thermal equilibrium one since in the Hamiltonian, the interaction part are not included. Since we assume three dimensional space is translational invariant, then the expectation value of the PNA depends on time x 0 and the initial time t 0 . As we will show later, the non-zero expectation value for the field φ 3 leads to the time dependent condensation which is the origin of the non-equilibrium time evolution of the system.
Next, one can write the expectation value of the PNA in Eq. (14) in terms of Green's function and the mean fields 1 ,φ [17], Eq. (17) can be understood by writing original field as a sum of mean field and quantum field, φ =φ + ϕ. The upper indices of the mean field and Green's function distinguish two different time paths in closed time path formalism [16]. Then Eq.(14) alters into 2 , where we have used Eq.(12) and the following relations, where τ is the Pauli matrix. 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.
A. 2PI Formalism in Curved Space-Time 2PI CTP EA in curved space-time has been investigated in [18] and their formulations can be applied to the present model. In 2PI formalism, one introduces non-local source term denoted as K and local source term denoted as J, where c ab is the metric of CTP formalism in [16] and c 11 = −c 22 = 1 and c 12 = c 21 = 0. The Legendre transformation of W leads to the 2PI EA [17], where S is the action written in terms of background fields, 3 and Γ Q is given by,

B. 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. We first derive SDEs for the field. The variation of the 2PI EA in Eq. (22) with respect to the scalar fieldφ leads to, and one obtains the following equation of motion of the scalar fieldφ, where the Laplacian in Friedman-Lemaître-Robertson-Walker metric is given by, The variation of the 2PI EA with respect to Green's function G leads to, Taking variation of Eq. (22) with respect to Green's function, one obtains, The second and third terms of above expression are computed using action in Eq.(23) and Γ Q in Eq.(24), respectively. Furthermore, Eqs.(28) and (29) lead to two differential equations, where x = ∇ µ x ∇ x µ and y = ∇ µ y ∇ y µ , and we have defined matrix multiplication as follows, 4 Next, we rescale Green's function, field and coupling constant of interaction as follows, where a t 0 stands for the initial value for the scale factor and we have defined a t 0 := a(t 0 ).
By using these new definitions, SDEs in Eqs.(26), (30) and (31) are written as, where we have defined, and have used the definition of shorthand notations which are given by, Note that the first derivative with respect to time which is originally presented in the expression of Laplacian, Eq.(27), is now absent in the expression of SDEs for the rescaled fields and Green's functions.

C. The initial condition for Green's function and field
In this subsection, the initial condition for Green's function and field and their time derivatives is determined. For this purpose, we consider the local source term J and the non-local source term K in Eq. (21). K is non-zero only if both x 0 and y 0 are equal to the initial time [16] and it has the following form [14], while J has the following form, j and κ are related to the matrix element of the density operator in Eq.(15) as follows, where C is a normalization factor which is determined so that the density operator satisfies the normalization condition Tr(ρ) = 1.
is the Euclidean action for classical path which corresponds to the Hamiltonian in Eq. (16). It is given by, where κ ab ij (−k) are given as [14], and . By comparing Eq.(45) and Eq.(46), one obtains j and κ(x) as follows, The second term of Eq.(47) is independent of field and it can be absorbed in the normalization factor. Note that the normalization factors in the numerator and denominator of Eq.(46) will be canceled out. Then, Eq.(46) alters into, whereS Now let us compute the initial condensation of the fieldφ(t 0 ) =φ(t 0 ) (see Eq.(32)). It is given as 5 , where we have computed the last term of Eq.(53) using Eq.(50). By applying Eq.(48) to Eq.(50), we can show the following relation, D(r) in Eq.(54) is defined as [14], To proceed the calculation, we denote J(x) as J(x) = 2a 3 t 0 j 1 . Then the initial condensation of field φ(t 0 , x) is given by, The initial condition for Green's functions for quantum fields is obtained by applying the method in [14]. The quantum field is defined so that its condensation vanishes by subtracting the condensation from the original field. In reference [14], the initial condition for Green's functions for the fields with vanishing condensation is derived. From Eq.(46), the density operator for the quantum fields φ a i − v i has the same form as the one in [14]. Therefore the initial condition for Green's functions for quantum fields is the same as the one obtained in [14]. We summarize the initial conditions for both Green's functions and condensations, Next we derive the time derivative of the field and Green's function at the initial time t 0 .
First we integrate the field equation in Eq.(35) with respect to time. By setting x 0 = t 0 , we obtain, Similarly, we integrate Eq.(36) with respect to time x 0 . By setting both x 0 and y 0 equal to t 0 , we obtain the following initial condition, Finally, we integrate Eq.(37) with respect to time y 0 . By setting both x 0 and y 0 equal to t 0 , we obtain another initial condition,

IV. 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. They are unambiguously written as a sum of two terms. The solution of the field is written as, ϕ i,free is the free part contribution whileφ i,int is the contribution due to cubic interaction.
Their expressions are given in Eqs.(A24) and (A25), respectively. The Green's function is also written as a sum of the free part and the interaction part, The expression of the free part is given in Eq.(A42) while the interaction partĜ ij,int is given in Eq.(A45).
In We first write the solutions of fields as, and where W i,k is defined as, f i,k and g i,k are the solutions which satisfy the following homogeneous differential equations, where Next we write down the solution of Green's function as follows, where ǫ ab is an anti-symmetric tensor and its non-zero components are given as ǫ 12 = 1 while θ(t) denotes a unit step function, where Q o(A) , R o(A) and E k are given as, where κ ab ij (k) is given in Eq.(48). corrections as, 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 around t 0 for 0 < t 0 x 0 as follows, Then one can keep only the following terms, and a (n) (x 0 ) for (n 2) are set to be zero. a (0) corresponds to the constant scale factor and a (1) (x 0 ) corresponds to linear Hubble parameter H(t 0 ). Thus it can be written as, where H(t 0 ) is given by, and t 0 > 0. Throughout this study, we only keep first order of H(t 0 ) as the first non-trivial approximation. For the case that Hubble parameter is positive, it corresponds to the case for the expanding universe. Under this situation,ȧ(x 0 ) = a(t 0 )H(t 0 ) andä(x 0 ) = 0. In these approximations, the second term of Eq.(39) is apparently vanished. Sinceȧ(x 0 ) is proportional to H(t 0 ), the third term of Eq.(39) involves second order of H(t 0 ). Hence, one can neglect it and the Riemann curvature R(x 0 ) is also vanished. Therefore,m 2 i (x 0 ) is simply written asm 2 i . Nowm 2 i are given as, Next we define ω i,k as, We consider Ω i,k (x 0 ) defined in Eq.(38). One can expand it around time t 0 as, where we have used Eqs. (66) and (75). Following the expression of the scale factor in Eq.(81), K is also divided into the part of the constant scale factor and the part which is proportional to H(t 0 ) as shown in Eqs.(C4) and (C5). In the above expression, H(t 0 ) is also included inÂ(t). Since we are interested in the PNA up to the first order of H(t 0 ), we expand it as and the derivative ↔ ∂˙acts on the first argument ofK and defined as follows, Each term of the PNA shown in Eqs. (91) and (92) can be understood as follows. The first term is the PNA with the constant scale factor. The second term with a prefactor is called the dilution effect. The third term with a is called the freezing interaction effect. The fourth term which corresponds to j 0 (x 0 ) 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 Eq.(89). It implies that the strength of the cubic interactionÂ(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 Eq.(38), as the scale factor grows, the physical wavelength becomes large. Therefore, the momentum and the energy of the particles becomes 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 factor H t 0 = 0.
Using Eq.(D6) with H t 0 = 0, one obtains it at the initial time where n i is the distribution functions for the Bose particles, Because we assumem 1 <m 2 , one obtains inequality n 2 < n 1 . From the expression above, the production rate of PNA at the initial time is negative for v 3 A 123 > 0. One also finds the rate is logarithmically divergent for the momentum (k) integration, where µ = O(m i ) (i = 1, 2) and k max is an ultraviolet cut off for the momentum integration.
With the expression, one expects that for the positive v 3 A 123 , the PNA becomes negative from zero just after the initial time and the behavior will be confirmed in the numerical simulation.

V. 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 constant A 123 and the initial value of the fieldφ 3,t 0 , 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 a t 0 , without loss of generality, one can set this dimensionless factor is as unity. For the other parameters of the model, we choosem 2 , B and ω 3,0 =m 3 as independent parameters since the massm 1 is written as, The temperature T and the expansion rate H(t 0 ) 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 Eq.(81), H(t 0 ) 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, .
The time dependence of PNA is plotted as a function of the dimensionless time defined as, where ω r 3,0 is a reference frequency. In terms of the dimensionless time, the condition of Eq.(98) is written as, How the PNA behaves with respect to time is discussed in the following subsection (V A-V C). 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 V D, we assign the unit for the parameters and estimate the ratio of the PNA over entropy density. the fixed time t. As the expansion rate becomes larger, the size of PNA becomes smaller.
B. The PNA with the shorter period Now we investigate the PNA with the shorter period. In Fig. 6, 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. The Fig. 7 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 differencẽ m 2 −m 1 is larger than the frequency ω 3,0 . Since B 2 , proportional to mass squared differencẽ 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 differencem 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. In Fig. 8, 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, 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 9 shows the dependence of expansion rate (H t 0 ). In this plot, the PNA gradually decreases as the expansion rate increases.

C. 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 Fig. 10, 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 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, The above equation is derived as follows. The Einstein's equations without cosmological where G is the Newton's constant. ρ is the energy density for radiation and it is given by, where ρ 0 is the initial energy density and we set a t 0 = 1. By setting x 0 = t 0 in Eq.(103), the initial Hubble parameter is given by, Then using Eq.(105), Eq.(103) becomes, Solving the equation above, one can obtain Eq.(102).
From the expression in Eq.(102), one needs to specify the unit of the Hubble parameter at t 0 . Through Eq.(105), 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 (t 0 ), one can write ρ 0 as follows, Hereafter, we assume that the temperature of the radiation T (t 0 ) 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, where M Pl is the Planck mass, M Pl = 1.2 × 10 19 (GeV). Then one can write the temperature T in GeV unit as follows, In the numerical simulation, the ratio H t 0 /T is given. Therefore, for the given ratio and g * , the temperature T in terms of GeV unit is determined. Then H t 0 in GeV unit also becomes, The masses of the scalar fieldsm i (i = 1, 2, 3) can be also expressed in GeV unit as, where we use the ratiosm i Ht 0 given in the numerical simulation.
As an example, we study the implication of the numerical simulation shown in Fig.10 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 Eq.(109). 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 H t 0 /T (t 0 ) = 10 −5 in Fig.10 to Eq.(109), one obtains T ∼ 10 13 (GeV) and H(t 0 ) ∼ 10 8 (GeV). The mass parameters are different Mass parameter (GeV) The shorter period The longer period m 1 2 × 10 11 4 × 10 9 m 2 3 × 10 11 5 × 10 9 ω 3,0 3.5 × 10 10 3.5 × 10 8 between the longer period case (the dotted line) and the shorter period case (the solid line).
One can also estimate the size of PNA. Here, we consider the maximum value of the PNA for the longer period case in Fig.10. We evaluate the ratio of the PNA over entropy density s, where s is given by, Substituting the temperature T = 10 13 (GeV) into Eq.(112), one obtains, From the equation above, we can achieve the ratio as 10 −10 by taking A 123 = 10 8 (GeV) and v 3 = 10 11 (GeV).

VI. 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 with U(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 an another source of particle number violation, the U(1) symmetry breaking mass term for the complex scalar is The increase of volume of the universe due to expansion, Freezing interaction The decrease of the strength of the cubic interactionÂ asÂ 123 − A 123 .

Redshift
The effective energy of particle as indicated in Eq.(38), introduced. The initial value for the condensation of the neutral scalar is non-zero. Using The results showed that the PNA depends on the interaction coupling A 123 and the initial value of the condensation of the neutral scalarφ 3,t 0 . It also depends on the mass squared difference. We found that the interaction coupling A 123 and the mass squared difference play a key role to give rise to non-vanishing PNA. Even if the initial value of the neutral scalar is non-zero, in the vanishing limit of interaction terms and the mass squared difference, the PNA will vanish. Another important finding was 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 rests which are the first order term contribute according to their origins. Those are summarized in Table III.
We have numerically calculated time evolution of the PNA and have investigated its dependence on the temperature, parameter B, 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 non-zero PNA is generated, it starts to oscillate. The has not changed. The magnitude also has not changed significantly during the same period.
If the period of the oscillation is long enough, the quarter of the period, from the origin to the first peak of oscillation, is also long enough. We are supposed to be around the peak at present. It also takes a long time for the asymmetry to be washed out during the next quarter of the period.
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 be g * ≃ O(100). Then when the initial temperature of the scalar fields is the same as that of the light particles, the simulation with Ht 0 T = 10 −5 corresponds to the case that the temperature of the universe is 10 13 (GeV) which is slightly lower than GUT scale ∼ 10 16 (GeV) [19,20]. The masses of the scalar fields in Fig.10 are different between the shorter period case and the longer period case as shown in Table II. In the shorter period case, the mass spectrum of the scalar ranges from 10 10 (GeV) to 10 11 (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 (A 123 ) and the initial expectation value (v 3 ).
Compared with the previous works [13,14,21,22], instead of assuming the non-zero 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 feature 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.
As the future work, the relation between the PNA and the observed BAU should be studied in radiation dominated era. Since nuclei are produced from the BAU, then the baryogenesis should be completed by the BBN. If the PNA is a seed of the BAU, then it should be created before the BBN. One also needs to consider the mechanism how the created PNA is transferred to the observed BAU. In this subsection, we provide the general solution of SDEs. Let us introduce the following differential equation for a field ϕ, where S is an arbitrary function of time x 0 and we will find the solution within the time range from x 0 = t 0 to x 0 = T . At the boundaries x 0 = t 0 and x 0 = T , we introduce the source terms of the form of delta function. The strength of the delta function is denoted as E(t 0 ) and F (T ), respectively. One may assume that field vanishes at x 0 < t 0 , One integrates Eq.(A1) with respect to time x 0 from t 0 − ǫ to t 0 + ǫ and obtains initial condition for the first derivative of the field as, where we have used the above assumption and taken limit ǫ → 0. t + 0 denotes t 0 + 0.
The method of variation of constants has been employed to determine the solution of Eq.(A1) and it is written as, where f (x 0 ) and g(x 0 ) are two linear independent solutions of the following homogeneous equation, Firstly, the following condition is imposed, Eq.(A4) becomes the solution of the differential equation Eq.(A1) when C i (x 0 ) satisfies another condition,Ċ The remaining task is to find C i (x 0 ) which satisfy the conditions in Eqs.(A6) and (A7). It is convenient to write these conditions in matrix form as, where we have defined, and it is a constant with respect to time. Integrating Eq.(A8) with respect to time from t + 0 to x 0 , one obtains, Therefore Eq.(A4) becomes, In this regard, it enables us to define, Then Eq.(A12) alters into To determine C i (t + 0 ), it is required the first derivative of field, Consequently, using definition Eq.(A13), the term which includes S(t 0 ) vanishes after integration. From now on, we will denoteφ(x 0 ) as ∂ϕ(x 0 ) ∂x 0 .
Let us now write C i (t + 0 ) in terms of ϕ(t 0 ) andφ(t 0 ). We consider the following equations, Next, one can write C i (t + 0 ) in matrix form as, Thus one substitutes these C i (t + 0 ) into Eq.(A14) and obtains, where we have used the initial condition in Eq.(A3) andK ′ [x 0 , y 0 ] is defined as, In the next subsection, we will use the obtained solution to provide the solution of SDEs.

The SDEs for the field
Next we move to consider the SDEs for the field in Eqs.(35). It is rewritten as, where we have defined S d i,x 0 as, To solve Eq.(A21), one first sets E(t 0 ) = 0 and F (T ) = 0 in Eq.(A1). Then, the general differential equation is similar to the SDEs for the field. The SDEs of the field in the form of integral equation is given by,

The SDEs for the Green's function
In this subsection, we consider the SDEs for Green's function in Eqs.(36) and (37). They are simply rewritten as, where we have defined, and E ac ik,k is given in Eq.(78). In the following, we will obtain SDEs for Green's functions at (x 0 , y 0 ) in the form of integral equation. Starting with the initial condition for Green's function at (t 0 , t 0 ), we obtain two expressions for the Green's function at (x 0 , y 0 ). The two expressions correspond to two paths shown in Figure 11 which are used to integrate the differential equation in Eqs.(A26) and (A27). They are given by 6 , G go x 0 t · Q ty 0 dt where the upper indices "br"and "go"denote the blue red path and the green orange path, respectively. Now let us explain how one can derive Eq.(A31). The steps are summarized below.
• We first consider the differential equation in Eq.(A26) which Green's functionĜ(t, t 0 ) 6 Dot multiplication describes matrix product corresponding their indices.
on the blue line (t 0 ≤ t ≤ x 0 ) satisfies. Using the solution of the general differential equation in Eq.(A19), we obtain the expression, whereĜ t 0 t 0 denotes the initial condition.
• Next we consider the differential equation in Eq.(A27) which Green's functionĜ(x 0 , t) on the red line (t 0 ≤ t ≤ y 0 ) satisfies. Using the solution of the general differential equation in Eq.(A19), we obtain the expression, whereĜ x 0 t 0 denotes the initial condition.
The Eq.(A32) is obtained through the steps similar to the above. The difference is as follows.
We first integrate the differential equation on the green line with the initial condition at (t 0 , t 0 ) and obtain the expression forĜ t 0 y 0 . Using it as the initial condition, we integrate the differential equation on the orange line and obtain the expression in Eq.(A32).

The derivation of free part for Green's function and its path independence
In the following subsection, we derive the free parts of Green's function which are the zeroth order of cubic interaction. From Eqs.(A31) and (A32), we can write them respectively as,Ĝ br 11. Two paths to obtainĜ(x 0 , y 0 ). We show the paths for the case x 0 < y 0 .
Both of the above expressions satisfy the differential equations in which we turn off the interaction part, namely, Below we show both expressions in Eqs.(A35) and (A36) lead to a single expression. Using Eqs.(59),(78) and (A28), we can rewrite them as follows, By using the following relation, we can show that two expressions are identical to each other and they can be summarized into a single expression. It is given by, where we have defined ǫ as, and θ(t) is defined as,

The interaction part of Green's function and its path independence
Now we move to consider the interaction parts of Green's function. From Eqs.(A31) and (A32), we can temporary define them respectively as, G go x 0 y 0 ,int := where Q and R are written in terms of the sameĜ andφ in Eqs.(A29) and (A30). Below, we will show the two expressions are the same to each other. Remind us that Q and R in Eqs.(A45) and (A46) are written in terms ofĜ x 0 y 0 ,int in Eq.(64) through the following differential equation, Substituting these expressions to Eqs.(A45) and (A46), we obtain, respectively. To show the equalities of Eqs.(A50) and (A52), we have used the following relations,Ĝ We complete the proof of equality of two expressions given in Eqs.(A45) and (A46). Since they are identical each other, from now on, we will useĜ br x 0 y 0 ,int in Eq.(A45). To summarize this subsection, let us write the SDEs of Green's functions in the form of integral equations. Omitting the upper and lower indices i, j, a and b, the interaction part of Green's function is written as, where Q and R are given by, Appendix B: Derivation for f (x 0 ) and g(x 0 ) up to first order of H(t 0 ) As was discussed in Appendix A, f (x 0 ) and g(x 0 ) are the solutions of a homogeneous differential equation. In this appendix, we derive those solutions for the case that the scale factor is given in Eq.(81). We present them within linear approximation with respect to omit the lower indices, i and k), The above equation leads to the following leading equations, As for the solution of Eq.(B11), we choose, Next one needs to solve f (1) (s). The solution is written in terms of linear combination of sine and cosine functions, where their coefficients C i depend on time. Since we can impose the following condition, one can show that C i (s) satisfy, where C ′ i (s) are defined as, From Eqs.(B15) and (B16), one can write C ′ i (s) as, With the initial conditions C i (s 0 ) = 0, C i (s) are written as, Now using Eqs.(B17) and (B18), the solution of f (1) (s) is written as, To summarize this part, let us write f (0) (s) and f (1) (s) in terms of original dimensional parameters. They are given by, Now we move to compute for g(x 0 ). In this regard, g(s) satisfies the same equation in Eq.(B10) which f (s) satisfies. The difference is that g 0 (s) is a cosine function, Applying the same procedure which we have used for the derivation of f 1 (s), we obtain, Finally, one rewrites g (0) (s) and g (1) (s) in terms of original variables. They are given by, In this appendix, we presentK i,x 0 y 0 ,k given in Eq.(68) within the linear approximation with respect to H(t 0 ). For simplicity, momentum index k is suppressed.K i,x 0 y 0 ,k is also expanded up to the first order with respect to H(t 0 ), namely, where we have defined, One can show that W is written in terms of the zeroth order solutions f   is written in terms of original parameters as, We defineK i andK ′ i as,K i [x 0 , y 0 ] := ThenK ′ i in Eq.(A20),K i andK ′ i are given in terms of original parameters by, i [x 0 , y 0 ] +K (1) i [x 0 , y 0 ], K Below we first carry out time integration of j 0 (x 0 ) 1st,A . One defines the new variable of integration as, Then one obtains that, where we have defined ω ± 12,k as, ω ± 12,k := ω 1,k ± ω 2,k . (D9) The next task is to integrate Eqs.(D6) and (D8) with respect to spatial momentum. Using those equations, Eq.(D1) leads to the following expression, (D14) We carry out the momentum integration of the above expressions numerically.
leads to the following expression,