format: text

text:
# Universality in the pair contact process with diffusion

## I Introduction

Out of equilibrium systems may display phase transitions analogous to those found in their equilibrium counterparts. These transitions are classified into distinct universality classes, characterized by a set of critical exponents hinr00 ; marr96 . Particular attention has been paid to one-dimensional systems with transitions from an active state into absorbing states, i.e., frozen configurations from which the system cannot escape. For these systems, so far, only two distinct universality classes have been firmly established: the directed percolation (DP) and the parity conserving (PC) universality class. While the former is ubiquitous and found in several systems with very different dynamical rules, the latter is only known to occur when extra symmetries are present hinr00 ; marr96 .

A model which has attracted quite some interest recently, because it may indeed belong to a novel universality class howa97  is the so-called pair contact process with diffusion (PCPD). Despite a rather intense activity in the past couple of years carl01 ; hinr01 ; odor00 ; noh01 ; dick02 ; park02 ; hinr02a ; kock03 ; odor03 , the understanding of the PCPD and its relation with other known models is still unsatisfactory. In the fermionic version of the model—the one studied here—each lattice site is either occupied by a single particle ($`A`$) or empty ($`0`$). The reactions are

$`\begin{array}{ccc}\{\begin{array}{ccc}AA0& & AAA\\ 0AA& & AAA\end{array}& \mathrm{with}\mathrm{rate}& \frac{(1p)(1d)}{2},\\ AA00& \mathrm{with}\mathrm{rate}& p(1d),\\ A00A& \mathrm{with}\mathrm{rate}& d\end{array}`$ (6)
with $`0<d<1`$, $`0<p<1`$. The analysis of the critical properties of the PCPD has shown to be much more difficult than all similar models analyzed so far and several scenarios have been proposed. First, a similarity of the exponents $`\beta /\nu _{}`$ and $`z=\nu _{}/\nu _{}`$ (where $`\nu _{}`$ and $`\nu _{}`$ are the correlation length exponents along the time and space directions and $`\beta `$ the order parameter exponent) with those of the PC class carl01  was reported, although in the PCPD there is no conservation of parity. It was later suggested that the PCPD could belong to a new universality class with exponents close to, but different from the PC values hinr01 . It has also been argued that there could be two universality classes odor00  at small and large diffusion rates ($`d`$), or continuously varying exponents noh01 , or that scaling may even be violated dick02 . More recently the option of a slow crossover to DP was also discussed hinr02a , although the general belief is that the PCPD belongs to a novel universality class noh01 ; dick02 ; park02 ; kock03 ; odor03 . Field theoretical methods, which have been successfully applied to other reaction–diffusion models card96  failed so far to clarify the critical properties of the PCPD howa97 ; muno02 .

In this paper we present some insights into the PCPD. We show that accurate numerical results from Monte Carlo (MC) simulations and density matrix renormalization group (DMRG) calculations convincingly demonstrate that for sufficiently long times and system lengths the exponents crossover towards the DP values. Corrections to scaling are, however, rather strong and only an accurate extrapolation of the effective critical exponents allows to identify the final asymptotic critical behavior. We give evidence that in the near-asymptotic region the model shows effective exponents close to the PC class values, which suggests that the critical behavior of the system is described by two competing fixed points.

## II Results on bulk particles and pair densities

Our MC simulations exploit a technique known as multispin coding newm99 . The basic idea is that in a simulation of 64 systems, each with $`L`$ sites, the occupation of site $`i`$ in the $`k`$th simulation is stored in the $`k`$th bit of 64-bit word $`A[i]`$. To perform the reaction $`AA0AAA`$ in all 64 systems at a randomly chosen site $`i`$ and its neighbors, one logical operation $`A[i+1]=A[i+1](A[i]A[i1])`$ suffices, where $``$ and $``$ are the logical operations OR and AND, respectively. Other reactions require only slightly more elaborate logical operations. A direct implementation along these lines might result in 64 simulations, each of which statistically correct, but strongly correlated to each other because the site selection is shared. To alleviate this correlation without sacrificing efficiency, we employ random bit patterns that decide which reaction will be attempted in which system. The strong point of multispin coding is its efficiency. As illustrated above, only a few logical operations (each usually carried out in a single clock cycle without delay) suffice to update a site in 64 systems simultaneously. For each combination of $`d`$ and $`p`$ we simulated 64 systems with $`L=\mathrm{100\; 000}`$ sites over $`3\times 10^6`$ Monte Carlo time units, in about 15 hours on a single–processor workstation. We also simulated $`16\times 64`$ systems with $`L=\mathrm{100\; 000}`$ sites over $`10^7`$ Monte Carlo time units on a parallel computer, for $`d=0.5`$ and $`p=\mathrm{0.152\; 45}`$, our estimate for the critical point.

A standard procedure to obtain critical exponents by means of Monte Carlo simulations is to study the decay of the particle density $`\rho (t)`$ as a function of time $`t`$, starting from a random configuration of particles. At the critical point one has $`\rho (t)t^{\beta /\nu _{}}`$. To monitor the decay it is convenient to define the effective exponent $`\delta _{\mathrm{eff}}\mathrm{ln}\rho (t)/\mathrm{ln}t`$. Typically $`\delta _{\mathrm{eff}}`$ is plotted as a function of $`1/t`$. At the critical point, in the limit $`t\mathrm{}`$, it approaches a finite value ($`\delta _{\mathrm{eff}}\beta /\nu _{}`$), while it deviates upwards or downwards with respect to this value in the inactive and active phases, respectively. This criterion allows to estimate the critical point location and the ratio $`\beta /\nu _{}`$ hinr00 . Some care, however, has to be taken when corrections to scaling are particularly strong, for example when $`\rho (t)t^{\beta /\nu _{}}\left(1+ct^\gamma \right)`$ with $`\gamma <1`$ and $`c`$ a constant. In this case $`\delta _{\mathrm{eff}}`$ plotted as a function of $`1/t`$ approaches $`\beta /\nu _{}`$ with an infinite slope. The ideal situation would be to plot $`\delta _{\mathrm{eff}}`$ as a function of $`1/t^\gamma `$ as in that case the approach to the asymptotic value would be linear. Further on in this paper we will give numerical evidence that the correction-to-scaling exponent $`\gamma `$ is close to $`\beta /\nu _{}`$ (which is much smaller than 1). A natural choice is therefore to plot $`\delta _{\mathrm{eff}}`$ as a function of the particle density $`\rho `$ instead of $`1/t`$, to avoid infinite slopes.

Besides the decay of the particle density $`\rho `$, we also consider the decay of the pair density $`\rho ^{}AA`$. Before presenting the results for the critical exponents we analyze the behavior of the ratio $`\rho /\rho ^{}`$ at the critical point. Such a quantity, which is shown in Fig. 1 for $`d=0.5`$, approaches a constant value ($`2.45`$) asymptotically for long times. The most important consequence of this fact is that in PCPD one can extract the critical exponent $`\beta /\nu _{}`$ both from the decay of the particle and pair densities ratio . Numerically the ratio is much better behaved than the individual densities $`\rho `$ and $`\rho ^{}`$, as fluctuations in the individual densities are highly correlated and largely cancel each other in the ratio.

The data of Fig. 1 also provide an estimate of the leading correction term in the asymptotic limit $`t\mathrm{}`$, which appears to be of the type $`\rho /\rho ^{}C(1+Dt^\gamma ^{})`$ with $`\gamma ^{}0.33`$, as shown by the linear approach to the asymptotic behavior of the data when plotted versus $`t^{0.33}`$. This exponent is much larger than its equivalent $`\gamma `$ in the particle and pair densities, which implies that the leading corrections for $`\rho `$ and $`\rho ^{}`$ cancel in the ratio $`\rho /\rho ^{}`$. We tested that a similar cancellation also occurs in the process $`A3A`$, $`2A0`$, which belongs to the PC universality class hinr00 : the leading correction in $`\rho `$ and $`\rho ^{}`$ scales as $`t^{0.6}`$, whereas the leading correction in the ratio $`\rho /\rho ^{}`$ scales as $`1/t`$.

Figure 2 shows a plot of $`\delta _{\mathrm{eff}}`$ at $`d=0.5`$ as a function of $`\rho `$ and $`\rho ^{}`$, calculated in the PCPD from the decay of particles and pairs, respectively. Three different values for $`p`$ have been plotted around the critical point which we estimate as $`p_c=\mathrm{0.152\; 45}(5)`$. For $`p>p_c`$ (inactive phase) and $`p<p_c`$ (active phase) $`\delta _{\mathrm{eff}}`$ rapidly veers up and down, as expected. At the critical point, $`\delta _{\mathrm{eff}}`$ approaches the $`y`$ axis with a finite slope, indicating that the leading correction to scaling is most likely described by an exponent roughly equal to $`\delta `$ itself. A parabolic fit through the data yields as a common estimate for the critical exponent $`\beta /\nu _{}=0.17`$. Similar calculations were repeated for other values of the diffusion coefficient $`d=0.05`$, $`0.1`$, $`0.2`$, and $`0.9`$ (the inset of Fig. 2 shows the case $`d=0.1`$), with the same results which we can summarize with the estimate $`\beta /\nu _{}=0.17(1)`$. This value is consistent with the DP class exponent $`\beta /\nu _{}=0.159`$.

An alternative way to calculate critical exponents would be to use a linear fit of the densities versus time on a double-logarithmic scale. However, such fits are hazardous when dealing with very slow convergence, as is the case here, and may lead to wrong estimates for the critical exponents. We illustrate this in Fig. 3 which shows the plot of the average particle density as a function of the time in a double-logarithmic scale for $`d=0.5`$ and $`d=0.9`$. In the former case a straight line (dashed) fits extremely well the critical density decay leading to the estimate $`\delta _{\mathrm{eff}}=0.219`$. The analysis of the effective exponent (see Fig. 2), however, provides a closer inspection of the local slopes of the double-logarithmic data. This analysis reveals some remaining curvature, and the final estimate of the exponent is significantly lower, compared to that obtained from the fit in the double-logarithmic scale. In the case of higher diffusion \[see Fig. 3(b)\] the curvature is more pronounced and clearly visible also in the double-logarithmic plot, which can be fitted by two straight lines with slopes $`\delta _{\mathrm{eff}}0.277`$ in the range $`4\mathrm{ln}t8`$, and with $`\delta _{\mathrm{eff}}0.212`$ for $`10\mathrm{ln}t15`$. Notice that the former exponent is consistent with that expected for the PC class ($`\delta _{\mathrm{PC}}=0.286`$ hinr00 ).

Again an extrapolation of the effective exponent (as done in Fig. 2) shows convergence to a value consistent with DP. Note that the value $`\beta /\nu _{}0.21`$ is consistent with the most recent Monte Carlo estimates for the PCPD odor02 ; kock03 ; hinr02a . In particular, Kockelkoren and Chaté, kock03  performed a series of Monte Carlo simulations for a bosonic version of the PCPD where the constraint of one particle per site is released. Their estimation of critical exponents is based on a straight-line fit to a double-logarithmic plot of $`\rho `$ versus $`t`$, from which they find $`\beta /\nu _{}=0.200(5)`$. This bosonic version of the model is claimed to suffer less from corrections to scaling than the fermionic case. Notice, however, that also in the fermionic PCPD studied here the density decay at $`d=0.5`$ (see Fig. 3(a)) is rather straight in a double-logarithmic plot for simulation times similar to those in Ref. kock03 . The advantage of the effective exponent analysis performed here is that it allows to extrapolate the numerical results to time scales beyond those actually simulated.

Next we present some DMRG results. DMRG karen  allows to calculate accurate stationary state probabilities for chains of moderate lengths carl99 . As usual in DMRG, we used open boundary conditions. In the PCPD on a lattice of finite length there are only two stationary states: a state with no particles and a state occupied by a single diffusing particle. To induce a finite density of particles we added a reaction $`0A`$ at the two boundary sites. The particle density decays from the two boundaries and forms a U-shaped profile. For chains of various lengths we calculated the density of particles $`\rho (L)`$ and of pairs $`\rho ^{}(L)`$ at the central site of a system of length $`L`$. At the critical point these quantities decay in the limit $`L\mathrm{}`$ as $`\rho (L)\rho ^{}(L)L^{\beta /\nu _{}}`$.

Figure 4 shows the effective exponent $`\gamma _{\mathrm{eff}}=\mathrm{ln}\rho (L)/\mathrm{ln}L`$ versus $`\rho `$ for $`d=0.5`$ at the critical point. As in Fig. 2 we include also the data for the pairs. Dotted and dashed lines are fits with polynomials in the densities. Again, a test of good convergence is that both exponents extrapolate to the same asymptotic value. This requirement seems indeed to be fulfilled and we find as extrapolation $`\beta /\nu _{}=0.27(4)`$. This exponent is again consistent with the DP value $`\beta /\nu _{}=0.252`$ hinr00 . Similar results have also been found for other values of the diffusion coefficient $`d`$. Extrapolations for $`d=0.2`$ are shown in the inset of Fig. 4. At small $`d`$ the maxima in $`\gamma _{\mathrm{eff}}`$ shift to longer $`L`$, thus extrapolations are somewhat less stable. In this case we take the estimate obtained from the particles $`\beta /\nu _{}=0.28(5)`$.

Previous DMRG results carl01  were restricted to the density of particles and to smaller systems than that studied here. The data for the effective exponent showed a monotonic behavior (except at very strong diffusion) and were analyzed using an extrapolation with polynomes in $`1/L`$. These extrapolations lead to a value consistent with the PC class exponent $`\beta /\nu _{}=0.50`$ carl01 . The present calculation, extended to the density of pairs and to longer systems, reveals that nonmonotonicity in the effective exponent is a common feature at all $`d`$. This nonmonotonicity leads to a rather strong decrease of the extrapolated exponent compared to the estimates of Ref. carl01 .

## III Results on the dynamical exponent

From the ratio of the exponents $`\beta /\nu _{}`$ and $`\beta /\nu _{}`$ one can estimate the dynamical exponent $`z=\nu _{}/\nu _{}`$. Since both $`\beta /\nu _{}`$ and $`\beta /\nu _{}`$ are consistent with DP, also the dynamical exponent $`z`$ agrees with the DP value $`z_{\mathrm{DP}}=1.58`$. It is however instructive to show the results of an independent calculation of $`z`$. This quantity can be obtained from a finite size scaling analysis of $`\mathrm{\Delta }`$, the gap of the Master operator, which is the inverse of the relaxation time of the system (see Ref. carl99  for details). As a function of the system length $`L`$ the gap decays as $`\mathrm{\Delta }L^z`$.

Figure 5 shows a plot of the effective exponent $`z_{\mathrm{eff}}=\mathrm{ln}\mathrm{\Delta }/\mathrm{ln}L`$ versus $`1/L`$. The calculations are similar to those reported in Ref. carl01 , but now for longer systems (up to $`L=46`$ compared to $`L=30`$ of Ref. carl01 ). The critical point locations were obtained from Monte Carlo simulations, which for this purpose are faster and more efficient than DMRG. Therefore we concentrated our computational efforts on a single value of $`p=p_c`$ and could obtain results for longer systems. As is clear from Fig. 5, the exponent $`z_{\mathrm{eff}}`$ is rather sensitive to the value of the diffusion rate $`d`$. As the estimates of $`\beta /\nu _{}`$ and $`\beta /\nu _{}`$ are instead rather stable as a function of $`d`$ we contribute this sensitivity to rather strong finite-size effects. Notice that the finite $`L`$ corrections change sign from the weak to the strong diffusion regime. The border value is around $`d=0.5`$ where $`z_{\mathrm{eff}}`$ has a very weak dependence on $`L`$. The data (see inset) run extremely close to the PC value $`z_{\mathrm{PC}}=1.75`$. At higher diffusion rates $`d0.80.9`$ the effective exponent $`z_{\mathrm{eff}}`$ for the range of sizes investigated is much lower than $`z_{\mathrm{PC}}`$. At the strongest diffusion investigated, extrapolations with different forms for the correction to scaling terms as $`1/L`$ or $`1/\sqrt{L}`$ yield values in the range $`1.5z1.65`$, which should be compared with the DP value $`z_{\mathrm{DP}}=1.58`$. Current Monte Carlo estimates from various authors hinr01 ; odor00 ; kock03  place the exponent $`z`$ in the range $`1.71.8`$ and the calculations were mostly performed in the weak diffusion regime.

We also performed a series of Monte Carlo simulations to calculate the exponent $`z`$ using finite-size scaling analysis. At the critical point and on a finite system the particle density decays as $`\rho =t^{\beta /\nu _{}}f(tL^z)`$, with $`f`$ a scaling function. For finite $`L`$, $`\rho `$ follows a power-law decay up to a characteristic time $`\tau `$ after which it drops exponentially. One expects that $`\tau `$ scales as $`\tau L^z`$. To estimate $`\tau `$ we calculated $`\rho `$ for lattices up to $`L=5000`$ and $`t=10^7`$ Monte Carlo time units. Figure 6 shows a plot of $`\mathrm{ln}(t^{\beta /\nu _{}}\rho )`$ versus $`\mathrm{ln}t`$ for various $`L`$ and $`d=0.9`$. The intersection of the data with a horizontal line at $`t^{\beta /\nu _{}}\rho =k`$ (with $`k`$ a constant) provides an estimate of $`\tau `$. As we work in a region where the particle density is rather low and fluctuations are large, and as the calculation of $`z`$ requires very smooth data, we performed averages over a large number of samples ($`>10^3`$). For the calculation we used $`\beta /\nu _{}=0.17`$, which is the value determined above, and $`k=2`$ (see Fig. 6). The inset shows a double-logarithmic plot of $`\tau `$ versus $`L`$ for $`d=0.2`$ and $`d=0.9`$ at their critical points. In the former case we restricted ourselves to $`L=2000`$ as the relevant times are typically longer at weak than at stronger diffusion, as expected. Notice that in both cases the data are well fitted by straight lines yielding the estimates $`z=1.70(3)`$ for $`d=0.2`$ and $`z=1.61(3)`$ for $`d=0.9`$, where the latter value is consistent with the dynamical exponent of directed percolation $`z_{\mathrm{DP}}=1.58`$. The results generally confirm the DMRG findings according to which the dynamical exponent is generically smaller, for finite $`L`$, at higher diffusivity. We also notice that by varying the value of $`\beta /\nu _{}`$ entering in the $`y`$ axis of Fig. 6 one changes the estimate for $`z`$. For instance, if we take $`\beta /\nu _{}=0.20`$, as calculated in Ref. kock03 , this leads to an increase of $`0.03`$ in the estimated value for $`z`$. The estimate of $`z`$ is rather stable for changes in the constant $`k`$.

## IV Results on surface densities

Boundary quantities are easily accessible in DMRG techniques carl97 , as one is basically forced to work with open boundary conditions. Surface criticality in absorbing phase transitions has been the subject of several studies in the past years both for models in the DP laur97  and in the PC laur98  universality classes. In the latter case, it is known that there are two distinct surface exponents depending on the type of boundary conditions applied laur98 . The results of a DMRG calculation of the surface critical exponents for a reaction-diffusion model in the PC class are presented in the Appendix. Here, we report on the surface critical exponent calculations for the PCPD, using the same types of boundary conditions as in the Appendix.

As in the calculation of the bulk particle density of the preceding section we inject particles through the reaction $`0A`$ at the boundary site labeled by the position $`i=1`$ in order to induce a finite density of particles in the system, and we measure the particle density $`\rho ^\mathrm{s}(L)`$ at the opposite boundary site $`i=L`$. Asymptotically for $`L\mathrm{}`$, we expect $`\rho ^\mathrm{s}(L)L^{\beta ^s/\nu _{}}`$, where $`\beta ^\mathrm{s}`$ is the order parameter surface exponent. The two different boundary conditions (BCs) applied at the site $`i=L`$ are: (a) No particles are allowed to leave the system from the boundary site and (b) particles may diffuse out of the system, i.e., the reaction $`A0`$ (with rate $`d`$) is added at that site. We refer to these as reflecting and absorbing boundary conditions, respectively.

In Fig. 7 we plot the effective exponent $`\gamma _{\mathrm{eff}}^s=\mathrm{ln}\rho ^\mathrm{s}(L)/\mathrm{ln}L`$ versus $`\rho ^s`$ in the case of reflecting BCs. Horizontal lines show the ratio $`\beta ^s/\nu _{}`$ for DP ($`=0.667`$ laur97 ; carl99 ) and PC ($`=0.72`$ laur98 ). In the DP case the different BCs produce the same critical exponent. Effective exponents in this case grow monotonically, contrary to what is found for bulk exponents. Notice that a cubic fit yields a quite stable estimate $`\beta ^s/\nu _{}=0.72(1)`$ in the range $`d0.5`$, a value actually consistent with the surface exponent for the PC class (see Appendix). Only at higher $`d`$ we observe some deviation from PC. The fact that the extrapolated surface exponents vary with $`d`$, while our current estimates for the bulk exponents are independent on $`d`$, is an indication that the former are not yet the correct asymptotic ones.

Figure 8 shows $`\gamma _{\mathrm{eff}}^s`$ versus $`\rho ^s`$ for the case of absorbing boundary conditions. Again for weak diffusion the exponent seems to extrapolate rather convincingly to values close to the PC class ($`\beta ^s/\nu _{}1.11`$, see Appendix), while for strong diffusion it increases to much larger values. Also in this case there is no clear signature of nonmonotonic behavior, except for the case $`d=0.5`$ where the data for the largest systems pass through a maximum.

We also analyzed the effective exponent data from the pair density $`\rho ^s`$ which are shown in the inset of Fig. 8 in the case of absorbing BCs, and plotted as functions of $`\rho ^s`$. In the range $`d0.5`$ the data extrapolate close to the PC surface exponent $`\beta ^s/\nu _{}1.11`$ as for the particle density. At very strong diffusion ($`d=0.9`$) the surface effective exponent shows a nonmonotonic behavior with a maximum around $`\gamma _{\mathrm{eff}}^s2.2`$. Notice that particle and pair exponents in this case are rather far apart from each other and it is quite hard to find a common extrapolation value. We would expect for $`\gamma _{\mathrm{eff}}^s`$ a similar behavior as for the bulk exponents, i.e., an increase followed by a decrease towards the asymptotic value. We suspect that in the present surface exponent calculation the decreasing side has barely been reached. So we tend to distrust the extrapolation as estimates of the genuine asymptotic behavior. They rather provide some insight on the preasymptotic region and actually point to a similarity with PC surface exponents at weak diffusion.

## V Discussion

To conclude, by combining Monte Carlo and DMRG calculations we analyzed the critical properties of the pair contact process with diffusion. This model has been the subject of increasing attention in recent years. Although the debate around it has not yet been settled, the main belief is that the PCPD belongs to a novel universality class which differs from the known DP and PC classes.

In our opinion, however, the most plausible scenario for the PCPD is that it ultimately falls into the DP universality class. The asymptotic behavior is, however, masked by rather strong finite size and time effects, characterized by small correction-to-scaling exponents, as our Monte Carlo simulations for the decays of the particle density $`\rho `$, the pair density $`\rho ^{}`$, and the ratio $`\rho /\rho ^{}`$ have demonstrated.

The exponents $`\beta /\nu _{}`$ and $`\beta /\nu _{}`$ extrapolated both from $`\rho `$ and $`\rho ^{}`$ appear to be stable as functions of the diffusion constant $`d`$ and actually consistent with the DP class values. The data show a nonmonotonic behavior both in time and system size, which in our opinion points to a crossover phenomenon between two competing types of critical behavior. The surface exponents, which we also investigated, turned out to be instead rather sensitive to the value of $`d`$; a sign, in our opinion, that the extrapolated values are probably not the true asymptotic ones. Interestingly enough, particularly at weak diffusion, the extrapolated values are rather stable and consistent with those expected for the PC class.

In early numerical studies of the PCPD carl01 ; hinr01 ; odor00 , restricted to shorter simulation times and system lengths compared to those considered here, several quantities as $`\beta /\nu _{}`$, $`\beta /\nu _{}`$, and $`\nu _{}/\nu _{}`$ were found to be quite consistent with the PC class values. It is now generally agreed that the PCPD does not belong to the PC universality class, as more extensive simulations performed by several groups have shown convincingly noh01 ; park02 ; dick02 ; hinr02a ; kock03 . Still one would like to understand if the observed similarity with the PC exponent is purely fortuitous or if there is some deeper reason for it. In our opinion the evidence given above that also the surface exponents extrapolate towards PC values in an intermediate regime strongly suggests that there is a genuine nonasymptotic PC-like regime, with a crossover to DP behavior at longer time scales.

A prototype model in the PC class is the branching annihilating random walk with even offsprings (BARWe), defined by the reactions $`A3A`$, $`2A0`$ plus diffusion hinr00 ; marr96 , which differs from the PCPD only for the reaction which creates particles. We argue that the early stages of the critical dynamics, when the system has a rather high particle density, are dominated by the annihilation process $`2A0`$, so that the substitution of the BARWe reaction $`A3A`$ with that of the PCPD $`2A3A`$ may result in a very weak perturbation of the system. Therefore a transient PC-like regime may be observed for $`t\tau _\mathrm{c}`$, where $`\tau _\mathrm{c}`$ is some crossover time. This argument may help to explain features observed in the PCPD, and should be equally valid for other models where the annihilation is of the type $`2A0`$ and with different creation rules $`nA(n+k)A`$ with $`n2`$, $`k>0`$; for such systems we expect a transient PC regime as well.

The study of reaction-diffusion systems where the annihilation and creation reactions involve $`n2`$ particles has recently drawn some attention kock03 ; park02b ; odor02 . In particular, we mention here the two cases recently considered by Ódor odor02  (i) $`3A5A`$, $`2A0`$ and (ii) $`4A5A`$, $`2A0`$. In model (i) he estimates $`\beta /\nu _{}0.28`$ (consistent with PC) for small diffusion rates and $`\beta /\nu _{}0.24`$ at stronger diffusion. Invoking some logarithmic corrections he claims that all values extrapolate to $`\beta /\nu _{}0.22`$ odor02 . In case (ii) the estimate is $`\beta /\nu _{}0.28`$ both at high and low diffusions odor02 , again consistent with the PC class value. The above observations suggest that these types of systems follow closely a critical behavior as described here for the PCPD, and it is thus plausible that they fall for sufficiently long times into the DP class. However, it may turn out to be quite difficult to show this numerically, as we expect that increasing the number of particles involved in the creation and annihilation reactions will lead to models even harder to simulate and analyze than the PCPD.

Very recently Kockelkoren and Chaté analyzed a similar set of models kock03 . In their formulation the fermionic constraint of only one particle per site is released. All the reactions of the type $`nA(n+k)A`$ and $`2A0`$ with $`n>2`$ were found to belong to the DP class. Surprisingly, in all those models the convergence to DP exponents seems to be quite fast (at least for $`\beta /\nu _{}`$) and not plagued by the strong corrections found in the fermionic models. It would be interesting to study the same models at different values of the diffusion constant, as in the PCPD the onset of crossover behavior is quite strongly influenced by the value of $`d`$.

We are grateful to J.D. Noh, H. Hinrichsen, M. den Nijs, and F. van Wijland for useful discussions.

*

## Appendix A Surface critical behavior in the parity conserving process 2A $``$ 0, A $``$ 3A

We present here some results on the surface critical behavior of the parity conserving process defined by the reactions $`2A0`$, $`A3A`$ and with single particle diffusion. We show how DMRG produces accurate surface critical exponents for this model, which are in good agreement with former Monte Carlo simulation results. For the single particle diffusion and pair annihilation we used the same rates as in Eq. (6), while we assign a rate $`(1p)(1d)`$ to the reaction $`0A0AAA`$. We restrict ourselves to a single value of the diffusion constant $`d=0.5`$.

We first estimated the critical point at $`p=p_c0.577(2)`$ by means of Monte Carlo simulations using a standard approach. As mentioned above, for surface universality in PC processes there are two possible types of boundary conditions leading to two distinct surface exponents laur98 . In the first case, the system is truncated at one edge and no particles are allowed to cross the boundary site; we refer to this as reflecting boundary conditions. In the second case, particles are allowed to drop from the boundary. We implemented this type of boundary condition adding the boundary reaction $`A0`$ (with rate $`d`$) which mimics the diffusion of particles out of the system. We refer to this implementation as absorbing boundary conditions.

Figure 9 shows the effective surface exponent $`\gamma _{\mathrm{eff}}^\mathrm{s}`$ versus $`1/L`$ in the case of reflecting boundary conditions, calculated both from the particle (circles) and pair (squares) densities. The same quantities are plotted in the inset in the case of absorbing boundary conditions. Notice that indeed the results confirm the existence of two distinct sets of surface exponents and that the data from pairs and particles merge for sufficiently long chains, indicating that both quantities decay with the same exponent. Our estimates $`\beta ^s/\nu _{}0.720(2)`$ in the former case and $`\beta ^s/\nu _{}1.10(1)`$ in the latter are obtained from a polynomial extrapolation in $`1/L`$. As the finite-size effects are rather small (see Fig. 9), the extrapolated values are not very sensitive to the type of correction to scaling term used in the extrapolation.

The Monte Carlo simulation results laur98  for the critical exponents are $`\beta ^\mathrm{s}=1.34(2)`$ and $`\beta ^\mathrm{s}=2.04(2)`$, for inactive and active boundary conditions, respectively. Combining these results with the PC class correlation length exponent $`\nu _{}=1.83(3)`$ hinr00 , one finds $`\beta ^s/\nu _{}=0.73(1)`$ (reflecting BCs) $`\beta ^s/\nu _{}=1.11(1)`$ (absorbing BCs), in very good agreement with the DMRG calculations.


text:
# What is the origin of chirality in the cholesteric phase of virus suspensions ?

Permanent address: \]Centre de Recherche Paul Pascal, UPR 8641, Avenue Albert Schweitzer, F-33600 Pessac, France.

## Abstract

We report a study of the cholesteric phase in monodisperse suspensions of the rod-like virus fd sterically stabilized with the polymer polyethylene glycol (PEG). After coating the virus with neutral polymers, the phase diagram and nematic order parameter of the fd-PEG system then become independent of ionic strength. Surprisingly, the fd-PEG suspensions not only continue to exhibit a cholesteric phase, which means that the grafted polymer does not screen all chiral interactions between rods, but paradoxically the cholesteric pitch of this sterically stabilized fd-PEG system varies with ionic strength. Furthermore, we observe that the cholesteric pitch decreases with increasing viral contour length, in contrast to theories which predict the opposite trend. Different models of the origin of chirality in colloidal liquid crystals are discussed.

preprint: HEP/123-qed

Chirality leads to a wonderful variety of liquid crystalline phases deGennes , from the cholesteric state to very complex structures such as smectic blue phases BPSm . In spite of its importance, the connection between chirality at the molecular scale and the macroscopic chiral structure of liquid crystalline phases is not yet understood, and quantitative prediction of the twist periodicity, also called cholesteric pitch, based on molecular features remains a unsolved problem. The elaboration of a rigorous statistical theory for the chiral interaction is one of the greatest challenges in the physics of liquid crystals.

We report and discuss experimental studies of the cholesteric phase of colloidal suspensions of filamentous bacteriophages, denoted fd. The rod-like virus fd is a micron-length semi-flexible polyelectrolyte formed by a single stranded DNA around which coat proteins are helicoidally wrapped. Therefore fd is a chiral object exhibiting a helical charge distribution. Suspensions of fd rods in aqueous solution form several liquid crystalline phases with increasing concentration, and especially a cholesteric or chiral nematic phase which lies between the isotropic and smectic phases PRLDogic . However, molecular chirality does not guaranty the existence of macroscopic chiral structure such as a cholesteric phase. For instance other chiral viruses such as Tobacco Mosaic Virus (TMV) TMV , or Pf1 langmuir  with a helical structure extremely similar to fd, exhibit only a nematic phase. Up to now, theory has failed to answer even the simple qualitative question: When does macroscopic chirality, such as twist, occur in liquid crystals formed by helical molecules? In this Letter, by studying a well characterized system of rod-like virus stabilized by grafted neutral polymer, we highlight that counterions play a crucial role. More specifically we show that ionic strength is one of the major physical parameters controlling the chirality of colloidal liquid crystalline suspensions.

The bacteriophage fd is a molecule with a contour length of L=880 nm, a diameter of $`D_{\text{bare}}`$=6.6 nm, a persistence length of $`\xi `$=2200 nm, a molecular weight of MW=1.64$`\times `$107 g/Mol, and a charge density of 10 e/nm in water at pH=8.1 langmuir . About 2700 copies of a single protein form the coat helicoidally wrapped following a 5-fold rotation axis combined with an 2-fold screw axis around the inner DNA makowski . The fd virus was grown and purified following standard biological protocols maniatis . To vary ionic strength (I), viruses were extensively dialyzed against a 20 mM TRIS-HCl buffer at pH=8.1 with an adjusted amount of NaCl. All measurements were done at room temperature and the concentrations were determined using spectrophotometry with absorption coefficient of 3.84 cm2/mg at 269 nm.

Although the fd virus is charged in solution, it is possible to prepare sterically stabilized colloids by irreversibly binding neutral polyethylene glycol (PEG) of different molecular weights (5,000 or 20,000 g/Mol) to the surface of the virus RoyalSoc . End functionalized PEG molecules (Shearwater corp.) have been used to covalently attach the water soluble polymer to the amino groups on the coat proteins of fd that are exposed to the solution. The refractive index increment, dn/dc, has been measured to estimate the degree of coverage of fd-PEG virus dn/dc , and 195$`\pm `$30 polymers PEG-20k (400$`\pm `$40 for PEG-5k) are attached to each virus. This result corresponds to nearly complete coverage of the virus by the polymer. Indeed, by considering grafted polymers as spheres of radius of gyration Rg, about $`\pi (1+D_{\text{bare}}/2R_g)\times L/2R_g`$ polymers can be close-packed around a virus of diameter $`D_{\text{bare}}`$ (Fig. 1). In the case of PEG-20K (Rg= 7 nm RgPEG ), this model predicts about 250 grafted polymers per virus, which is close to measured value.

The isotropic - cholesteric (I-N\*) transition has been investigated for this system of fd virus grafted with PEG RoyalSoc . According to Onsager Onsager , an isotropic to nematic phase transition occurs in suspensions of hard rods when the rod number density $`n`$ reaches $`n\pi L^2D16`$ (eq. 1). This is experimentally observed for bare fd (without any grafted polymer) over a wide range of ionic strengths as shown with square symbols in Fig. 1. Note that there is no difference in Onsager’s theory between the isotropic-nematic and isotropic-cholesteric phase transitions. The free energy difference between the nematic and cholesteric phase is negligible compared to the free energy difference between the isotropic and nematic phase Saupe .

Because the fd virus is charged, it is necessary to account for electrostatic repulsion by considering the colloids as having an effective diameter ($`D_{\text{eff}}`$) larger than their bare diameter ($`D_{\text{bare}}`$) with $`D_{\text{eff}}`$ roughly proportional to the Debye screening length Onsager . $`D_{\text{eff}}`$ therefore decreases with increasing ionic strength Tang95  and approaches $`D_{\text{bare}}`$ at high ionic strength. By substituting $`D_{\text{bare}}`$ with $`D_{\text{eff}}`$ in equation (1), it follows that the concentration of fd at the I-N\* transition is inversely proportional to $`D_{\text{eff}}`$. For fd-PEG rods, there is an ionic strength at which the coexistence concentrations of the I-N\* transition become independent of ionic strength. The crossover from the regime where the rods are electrostatically stabilized to where they are sterically stabilized by “soft” repulsion between virus-bound polymers is shown in Fig. 1. This corresponds to a ionic strength of I= 3 mM for fd-PEG-20k and I= 20 mM for fd-PEG-5k. At this ionic strength, the diameter associated with electrostatic interactions ($`D_{\text{eff}}`$) equals the diameter of the virus plus grafted polymer. It is therefore possible to extract an effective diameter for the polymer coated virus. In this range of ionic strengths, between 15 and 200 mM, the polymer is in a good solvent and the effects of NaCl on PEG properties are negligible in terms of radius of gyration and viscosity bailey . From Fig. 1, the effective diameter of fd-PEG-20k at I= 3 mM is 37 nm, which is approximately equal to the hard diameter of fd ($`D_{\text{bare}}`$= 6.6 nm) plus twice the diameter of PEG ($`R_g`$= 7 nm RgPEG ): $`D_{\text{bare}}+4R_g=35`$ nm. The fd-PEG-5k complex has $`D_{\text{eff}}=17`$ nm at high ionic strength, while $`D_{\text{bare}}+4R_g=19`$ nm ($`R_g`$= 3 nm for PEG-5k RgPEG ).

The good agreement between the measured and estimated value of the fd-PEG effective diameter supports the model of polymer behaving mainly as a sphere attached to the surface of the virus. Thus, the independence of the phase diagram on ionic strength reported in Fig. 1 shows that fd virus coated with PEG-20k is a system whose stability is entirely determined by polymer repulsion. This result has been confirmed by measuring the nematic order parameter of fd-PEG-20k close to the I-N\* transition. Fig. 2 indicates that the order parameter, determined by birefringence under magnetic field, is independent of ionic strength just as is the phase diagram RXkirstin . Therefore, fd-PEG-20k liquid crystal turns out to be a sterically stabilized system, that, surprisingly, continues to exhibit a cholesteric phase as shown in Fig. 3. Indeed, one could expect that grafted polymers on virus surface would screen chiral interactions, leading thus only to a nematic phase. To our knowledge, fd virus coated with PEG is the first completely sterically stabilized system which still has chiral features.

The insert of Fig. 3 presents a typical “fingerprint” texture of fd-PEG suspensions characteristic of a cholesteric phase observed by optical microscopy between crossed polarizers. After a few days of equilibration following the preparation of the samples, we have measured the dependence of cholesteric pitch on concentration and ionic strength (Fig. 3). Paradoxically, fd-PEG-20k liquid crystals have a cholesteric pitch which depends on ionic strength. Close to the I-N\* transition where both the phase diagram and order parameter are independent of ionic strength (Figs. 1 and 2), the pitch varies by a factor 2 from 100 $`\mu `$m at I=15 mM to 200 $`\mu `$m at I=110 mM. In this range of low rod concentrations the mean interparticle distance ($`d_{\text{inter}}`$) is about 60 nm: the fd-PEG surfaces are then separated by typically 3 bare diameters RXkirstin . The corresponding volume fraction at the I-N\* transition is $`\mathrm{\Phi }`$=0.2 where $`\mathrm{\Phi }`$ is defined by $`\mathrm{\Phi }=cN{}_{A}{}^{}L\pi (D_{\text{eff}}/2){}_{}{}^{2}/M_W`$ with $`N_A`$ being the Avogadro’s number. We can then calculate the concentration (c\*) for which the grafted polymers of fd-PEG-20k begin to overlap (corresponding to $`\mathrm{\Phi }`$\*=1): this gives c\*$``$35 mg/ml, which is close to the concentration c0=28 mg/ml where the pitch is independent of ionic strength (Fig. 3). The concentration c\* separates two regimes: we mainly focus on the dilute regime (c$`<`$c\*), where chiral rods do not overlap. Complexity increases for c$`>`$c\* because of short-range interactions between the polymers as well as electrostatic interactions between the viruses, which can induce intricate behaviors similar to those seen for instance in DNA condensation PhysicsToday .

What then is the mechanism of transmission of chiral interactions between sterically stabilized viruses? The first model of the origin of intermolecular twist was proposed by Straley who considered two screw-like molecules with excluded-volume interactions Straley . He predicted the pitch to be independent of the concentration, which is not observed experimentally. Moreover, if chirality is transmitted by some microscopic detail of the virus, we then expect that the chiral effects will be diminished for fd grafted with polymers. However, the same range of pitch has been seen for the fd-PEG system as for bare fd (between 10 and 200 $`\mu `$m) langmuir . That seems to exclude the strict steric model of packing of screws, because the thread of the screws of fd virus is covered by a 30 nm thick polymer coat. Even with the assumption that the polymers grafted to the helical protein coat replicate the chiral symmetry of the underlying bare virus, it does not explain why the cholesteric pitch should depend on ionic strength.

If the chirality is not transmitted at the microscopic length scale of the helical periodicity of the protein coat (with an axial repeat of 3.3 nm for fd virus makowski ), the existence of a superhelical twist, where chirality occurs at a much larger length scale, may be proposed. Brownian motion causes the virus to bend back and forth, and due to the coupling between flexibility and the underlying structural chirality imposed by the helical arrangement of coat proteins, the virus could twist more in one direction than in another leading to the creation of a “cork-screw” shape with a helical pitch of order a micron. A test \[suggested by F. Amblard\] of this model would be to measure the distribution of angles between the tangent vectors of the ends of a single virus confined to a surface. This is a Gaussian distribution for the worm-like chain model, but will have a maximum at a non-zero angle for a bent molecule, such as a cork-screw of length less than one full pitch. However, the “cork-screw” model can not be the whole story as Pf1, a chiral and flexible virus with a structure very similar to fd, has only a nematic phase, which implies that not all chiral flexible molecules form cork-screws. Additionally, completely rigid rods like cellulose Revol  and chitin Revol2  form a cholesteric phase implying that the cork-screw distortion is not necessary for cholesteric production. An interesting experiment would be to coat these rigid molecules with long polymer. If they still form the cholesteric phase then the cork-screw model would be rejected.

One of the most advanced theories of cholesteric assemblies of charged rods has been proposed by Kornyshev et al. for DNA. They suggest that various details of the helical symmetry of macromolecules surface charge pattern may be responsible for macroscopic twist leikin . The main idea is based on direct chiral electrostatic interaction between long helical molecules. Their attempt is promising, but how can such a theory account, even qualitatively, for our results on fd virus stabilized with polymer? Moreover, the prediction of a cholesteric pitch increasing with the square of the viral contour length ($`PL^2`$) shown both by Kornyshev et al. leikin  and by Harris et al. Harris  is in contradiction to our results obtained with mutant viruses (without any polymer grafted). Indeed, we are able to grow monodisperse mutants RoyalSoc  where only the contour length varies, keeping unchanged the local structure of the viruses mutants . We find that the cholesteric pitch decreases with increasing viral contour length as shown in Fig. 5. These experiments suggest that current theories do not capture all the essential elements of the connection between microscopic and macroscopic chirality in colloidal liquid crystals.

In conclusion, we have presented in this Letter a system of chiral colloidal rods coated with neutral polymer, which behaves as a sterically stabilized suspension because the phase diagram and order parameter are independent of ionic strength. However, this liquid crystal has the baffling feature that it exhibits a cholesteric phase whose pitch depends on ionic strength. If the cholesteric phase is induced by the charges on the the virus which are buried by a polymer coat of thickness corresponding to several Debye screening lengths, then we conclude that the cholesteric phase arises from very weak and subtle forces. Our result highlights the fundamental role of counterions in the existence of a macroscopic chirality and should stimulate new theoretical investigations on this fundamental topic of the origin of chirality in liquid crystalline phases.

This research was supported by NSF/DMR grants. We acknowledge Kirstin Purdy for her assistance in the preparation of virus mutants.


text:
# Affleck-Dine mechanism with negative thermal logarithmic potential

## I Introduction

Baryogenesis is one of the main issues in the theories of the early universe. There have been proposed various mechanisms, among which Affleck-Dine (AD) mechanism AD  is a promising candidate in theories with supersymmetry (SUSY). One of the features which distinguish supersymmetric theories from ordinary ones is the existence of flat directions, of which the AD mechanism makes use. In fact, there are many flat directions in the minimal supersymmetric standard model (MSSM), along which there are no classical potentials. Such a flat direction is described by a complex scalar field $`\mathrm{\Phi }`$, called AD field. It is crucial to determine the dynamics of the AD field to estimate the resultant baryon asymmetry. However, it has been found that the AD baryogenesis is complicated by the two causes. One is the thermal effect thermal ; thermal2 ; anisimov , and the other is the existence of Q balls Kusenko ; Enqvist ; Kasuya1 .

As for the thermal effects, we must take account of those which appear at both one-loop and two-loop orders. The former contribution to the effective potential of the AD field is the thermal mass term which might cause the early oscillation of the AD field. The latter effect, on which we would like to focus in this paper, induces the logarithmic potential (which we call thermal log potential hereafter) and hence totally changes the dynamics if its coefficient is negative. In that case, the AD field might be trapped far away from the origin for a long time until the temperature of the universe becomes low enough. In particular, such a trap is likely to occur in the gauge-mediated supersymmetry breaking (GMSB) models gmsb  where the potential for the AD field becomes flat at large amplitude. Therefore, we assume the GMSB models in this paper. Due to the enduring trap by the negative thermal log potential, the dynamics of the system is so involved that little has been known about whether the AD baryogenesis really works. Although the sign of the thermal contribution to the potential at two-loop order depends on the choice of flat directions, it is expected to be negative for many flat directions. Hence it is important to consider the case when the negative thermal log potential is effective.

The other element which we need to incorporate is the existence of Q balls Kusenko ; Enqvist ; Kasuya1 . After the AD field starts to oscillate, it experiences spatial instabilities and deforms into nontopological solitons, Q balls. Especially, the Q-ball formation is inevitable for a flat potential as in the GMSB models Kasuya1 ; Kasuya2 ; Kasuya3 . Once Q balls are formed, they absorb almost all the baryon numbers Kasuya1 ; Kasuya3 . In the usual picture, the Q-ball formation takes place after the baryon number is fixed. However, as shown later, the Q-ball formation plays an essential role to fix the baryon number if the AD field is trapped by the negative thermal log potential.

In this paper, therefore, we consider the AD baryogenesis with the negative thermal log potential in the GMSB models, taking into account the Q-ball formation. Since the negative thermal log potential retards the oscillation of the AD field considerably, the dynamics is quite different from that for the ordinary AD baryogenesis, and we need to examine the scenario carefully. It is found that the existence of the Hubble-induced A-term spoils the baryogenesis through the Q-ball formation, while it is possible to explain the baryon asymmetry of the present universe in the absence of the Hubble-induced A-term.

## II Affleck-Dine Mechanism and Q-ball formation

In this section we briefly review the Affleck-Dine mechanism and properties of the Q balls. In the MSSM, there exist flat directions, along which there are no classical potentials in the supersymmetric limit. Since flat directions consist of squarks and/or sleptons, they carry baryon and/or lepton numbers, and can be identified as the AD field. These flat directions are lifted by both the supersymmetry breaking effects and the nonrenormalizable operators with some cutoff scale. In the gauge-mediated SUSY breaking model, the potential of a flat direction is parabolic at the origin, and almost flat beyond the messenger scale Kusenko ; Kasuya3 ; Gouvea ,

$$V_{gauge}\{\begin{array}{cc}m_\varphi ^2|\mathrm{\Phi }|^2\hfill & (|\mathrm{\Phi }|M_S),\hfill \\ M_F^4\left(\mathrm{log}\frac{|\mathrm{\Phi }|^2}{M_S^2}\right)^2\hfill & (|\mathrm{\Phi }|M_S),\hfill \end{array}$$
(1)
where $`m_\varphi `$ is a soft breaking mass $``$ O(1 TeV), $`M_F`$ the SUSY breaking scale, and $`M_S`$ the messenger mass scale. Since the gravity always exists, flat directions are also lifted by the gravity-mediated SUSY breaking effects EnqvistMcDonald98 ,

$$V_{grav}m_{3/2}^2\left[1+K\mathrm{log}\left(\frac{|\mathrm{\Phi }|^2}{M_G^2}\right)\right]|\mathrm{\Phi }|^2,$$
(2)
where $`K`$ is the numerical coefficient of the one-loop corrections and $`M_G`$ is the gravitational scale ($`2.4\times 10^{18}`$ GeV). This term can be dominant only at high energy scales because of small gravitino mass $`O(1\text{ GeV})`$.

The nonrenormalizable operators, in addition, lift flat directions. Assuming that there exist all the nonrenormalizable operators consistent with both the gauge symmetries of the standard model and the R-parity in the superpotential, most of flat directions are lifted by the following operator:

$$W=\frac{1}{nM^{n3}}\mathrm{\Phi }^n,$$
(3)
where $`M`$ is the cutoff scale. Not only does it lift the flat potential as

$$V_{NR}=\frac{|\mathrm{\Phi }|^{2n2}}{M^{2n6}},$$
(4)
but it supplies the baryon number violating A-term as

$$V_A=a_m\frac{m_{3/2}}{nM^{n3}}\mathrm{\Phi }^n+\mathrm{H}.\mathrm{c}.,$$
(5)
where $`a_m`$ is a complex constant of order unity, and we assume the vanishing cosmological constant. The schema of the zero-temperature potential $`V_{gauge}+V_{grav}+V_{NR}`$ is shown in Fig. 1.

For the AD baryogenesis to work successfully, it is necessary that the AD field has a large expectation value during the inflationary epoch 111Actually, the AD field only has to have a large expectation value before it starts to oscillate. However we assumed that it has a large expectation value during inflation for definite argument.. To this end we require the existence of the four-point coupling of the AD field to the inflaton $`I`$ in the Kähler potential, which leads to the negative Hubble-induced mass term:

$$V_H=c_HH^2|\mathrm{\Phi }|^2,$$
(6)
where $`c_H`$ is a positive constant of order unity. Also there might be a Hubble-induced A-term, which comes from the three-point coupling of the AD field to the inflaton in the Kähler potential. We write it as

$$V_{AH}=a_H\frac{H}{nM^{n3}}\mathrm{\Phi }^n+\mathrm{H}.\mathrm{c}.,$$
(7)
where $`a_H`$ is a complex constant of order unity. Though the presence of this term usually does not change the resultant baryon asymmetry, it does change the dynamics of the AD field when the negative thermal log potential traps the AD field, as shown later. Since the Hubble-induced A-term is not necessarily present, we will consider both cases with and without this term.

Let us consider the thermal effects on the effective potential for the AD field. During the oscillation of the inflaton, there is dilute plasma with temperature $`T(T_{RH}^2HM_G)^{1/4}`$, where $`T_{RH}`$ is the reheating temperature. If the AD field directly couples with fields $`\psi _k`$ in the thermal bath, it acquires a thermal mass term in the effective potential at one-loop order:

$`V_T^{(1)}`$ $``$ $`{\displaystyle \underset{f_k|\mathrm{\Phi }|<T}{}}c_kf_k^2T^2|\mathrm{\Phi }|^2,`$ (8)
$`=`$ $`f_{(1)}^2T^2|\mathrm{\Phi }|^2\mathrm{for}|\mathrm{\Phi }|<f_{(1)}^1T,`$
where $`c_k`$ is a real positive constant of order unity, and $`f_k`$ is a Yukawa, or gauge coupling constant between the AD field and $`\psi _k`$ thermal2 ; afhy . For simplicity, we omit the subscript $`k`$ hereafter with the understanding that $`f_{(1)}`$ denotes either Yukawa or gauge coupling constant. Note that this effect is exponentially suppressed by the Boltzmann factor when the effective mass of the thermal particle, $`f_{(1)}|\mathrm{\Phi }|`$, is larger than the temperature. In addition to this term, there is another thermal effect on the potential, which appears at two-loop order, as pointed out in Ref. anisimov . This comes from the fact that the running of the coupling constant is modified by integrating out heavy particles which directly couples with the AD field. This contribution is given by afhy2

$$V_T^{(2)}a_Tf_{(2)}^4T^4\mathrm{log}\frac{|\mathrm{\Phi }|^2}{T^2},$$
(9)
where $`a_T`$ is an O(1) constant, and $`f_{(2)}^4=g^4,y^4,y^2g^2`$ ($`g`$ and $`y`$ are the gauge and Yukawa coupling constants). The sign of $`a_T`$ depends on the flat directions, and it is expected to be negative for many flat directions, for instance, $`\overline{u}\overline{d}\overline{d}`$ direction, though it is complicated to determine the precise value of $`a_T`$. Without specifying flat directions, we would like to focus on the case with the negative $`a_T`$ in the following discussion. (The AD baryogenesis for the positive $`a_T`$ was considered in Ref. Kasuya3 .)

When the AD field starts coherent rotation in the potential, the baryon number is generated and fixed immediately in the ordinary AD mechanism, since the A-terms becomes ineffective soon after the oscillation of the AD field. Then the baryon number density is estimated as

$$n_B(t_{osc})\epsilon \omega \varphi _{osc}^2,$$
(10)
where $`\epsilon (1)`$ is the ellipticity parameter, which represents the strongness of the A-term, and $`\omega `$ and $`\varphi _{osc}`$ are the angular velocity and the amplitude of the AD field respectively, at the beginning of the oscillation (rotation) in its effective potential. However, this is not the case when the AD field is trapped by the negative thermal log potential given by Eq. (9), which will be discussed in the next section. In fact, the baryon number is not fixed even after the AD field starts coherent oscillation, if the Hubble parameter is much smaller than the mass scale of the radial direction. It should be emphasized that Eq. (10) is valid as long as the baryon number is fixed soon after the oscillation commenced.

Lastly, we comment on the Q balls. In the course of the coherent oscillation, the AD field experiences spatial instabilities, and deforms into nontopological solitons, Q balls Kusenko ; Enqvist ; Kasuya1 . When the effective potential is dominated by the zero-temperature term $`V_{gauge}`$, the gauge-mediation type Q balls are formed, whose properties are as follows Dvali :

$$M_QM_FQ^{3/4},R_QM_F^1Q^{1/4},$$
(11)
where $`M_Q`$ and $`R_Q`$ are the mass and the size of the Q ball, respectively. If the mass per unit charge, $`M_FQ^{1/4}`$, is smaller than the proton mass $``$ 1GeV, the Q ball is stable against the decays into nucleons, from which it follows that Q balls with very large $`Q`$ can be stable. From numerical calculations Kasuya1 ; Kasuya3 , it is known that Q balls absorb almost all the baryon charges which the AD field obtains, and the typical charge is estimated as Kasuya3

$$Q\beta \left(\frac{\varphi _{osc}}{M_F}\right)^4,$$
(12)
where $`\beta 6\times 10^4`$. Consequently, the present baryon asymmetry should be explained by the charges which come out of the Q balls through the evaporation and decay of Q balls.

In the case of the unstable Q balls, they decay into nucleons. The decay rate is given by Coleman

$$\left|\frac{dQ}{dt}\right|\frac{\omega ^3A}{192\pi ^2},$$
(13)
where $`A`$ is a surface area of the Q ball.

In the case of the stable Q balls, the evaporation is the only way to extract the baryon charges from Q balls. The total evaporated charge from the Q ball is estimated as Kasuya3 ; Laine ; Banerjee ,

$$\mathrm{\Delta }Q10^{15}\left(\frac{m_\varphi }{\text{TeV}}\right)^{2/3}\left(\frac{M_F}{10^6\text{GeV}}\right)^{1/3}Q^{1/12}.$$
(14)
Hence the baryon number density is suppressed by the factor $`\mathrm{\Delta }Q/Q`$, in comparison with the case of unstable Q-ball production.

## III Dynamics of Affleck-Dine field

In this section, we present a detailed discussion on the dynamics of the AD field. First we derive the condition that the AD field is trapped by the negative thermal log potential. In the following argument we set $`|a_m|=c_H=|a_H|=c_k=a_T=1`$ for simplicity. During inflation the Hubble parameter is roughly constant, and the initial amplitude of AD field is set far away from the origin due to the negative Hubble-induced mass term Eq. (6). The position of the global minimum is

$$\varphi (HM^{n3})^{1/(n2)},$$
(15)
where $`\varphi \sqrt{2}|\mathrm{\Phi }|`$. After inflation ends, the Hubble parameter decreases, and so does the minimum of the potential. Since the AD field tracks the instantaneous minimum, Eq.(15) still holds after inflation. The AD field is trapped by the negative thermal log potential when $`|V_T^{(2)}|>|V_H|`$ is satisfied. Here we assumed that the trap due to the negative thermal log potential occurs while the AD field tracks the instantaneous minimum given as Eq.(15). The temperature $`T_c`$ and the amplitude $`\varphi _c`$ at the beginning of the trap are

$`T_c`$ $`=`$ $`\{\begin{array}{cc}f_{(2)}^{\frac{n2}{n}}T_{RH}^{\frac{n1}{n}}M_G^{\frac{n1}{2n}}M^{\frac{(n3)}{2n}}\hfill & \mathrm{for}T_{RH}<T_{},\hfill \\ & \\ T_{}\hfill & \mathrm{for}T_{RH}>T_{},\hfill \end{array}`$ (17)
$`\varphi _c`$ $`=`$ $`\left(f_{(2)}^2T_c^2M^{n3}\right)^{\frac{1}{n1}},`$ (18)
where we defined

$$T_{}f_{(2)}^{n2}M_G^{\frac{n1}{2}}M^{\frac{n3}{2}}.$$
(19)
It should be noted that $`T_{RH}>T_{}`$ requires a very high reheating temperature which might be constrained by the gravitino problem Gouvea . However, if the Q balls dominate the universe, they must decay and reheat the universe with the decay temperature, $`T_dO(1\mathrm{G}\mathrm{e}\mathrm{V})`$. Thus the gravitino problem is alleviated due to the entropy production of the decay of Q balls in this case, and such a high reheating temperature might be possible. If it were not for the gravitino problem, the reheating temperature could be as large as $`T_{inf}10^{16}`$ GeV. We take this value, because the COBE data implies that $`V_{inf}^{1/4}/ϵ_{sr}^{1/4}6.7\times 10^{16}`$ GeV LyRi , where $`ϵ_{sr}<1`$ is the slow-roll parameter, and the instantaneous reheating is assumed for conservative discussion.

Next we consider when the trap due to the negative thermal log potential ends. The AD field starts an oscillation when $`|V_T^{(2)}|<V_{gauge}`$ is first satisfied, hence the temperature $`T_{osc}`$ and the amplitude $`\varphi _{osc}`$ at that moment are

$`T_{osc}`$ $`=`$ $`f_{(2)}^1M_F,`$
$`\varphi _{osc}`$ $`=`$ $`\mathrm{min}[\varphi _\alpha ,\varphi _\beta ],`$ (20)
where $`\varphi _\alpha \left(M_F^2M^{n3}\right)^{\frac{1}{n1}}`$ and $`\varphi _\beta M_F^2/m_{3/2}`$. Note that the amplitude of the AD field, $`\varphi _{osc}`$, at the beginning of oscillation is different from that in the usual scenario. In the ordinary AD mechanism, the AD field starts oscillating at $`\varphi _\alpha `$ or $`\varphi _\gamma (m_{3/2}M^{n3})^{\frac{1}{n2}}`$ depending on the form of the potential shown in Fig. 1, if one does not take account of the thermal effects. Such modification is due to the different dependence of $`V_H`$ and $`V_T^{(2)}`$ on $`\mathrm{\Phi }`$.

From the above discussion we can derive the condition that the AD field is trapped by the negative thermal log potential. It requires the following inequalities to be satisfied.

$`\varphi _c`$ $`>`$ $`\mathrm{Max}[\varphi _\alpha ,\varphi _\gamma ],`$ (21)
$`f_{(1)}\varphi _c`$ $`>`$ $`T_c,`$ (22)
$`\varphi _{osc}`$ $`>`$ $`M_S.`$ (23)
Then we arrive at the following conditions.

$`10^5\mathrm{GeV}`$ $`m_{3/2}`$ $`m_{3/2}^{(u)},`$
$`10\mathrm{TeV}`$ $`M_F`$ $`M_F^{(u)},`$
$`T_{RH}^{(l)}`$ $`T_{RH}`$ $`T_{RH}^{(u)},`$ (24)
where $`m_{3/2}^{(u)}`$, $`M_F^{(u)}`$, $`T_{RH}^{(l)}`$ and $`T_{RH}^{(u)}`$ are given as

$`m_{3/2}^{(u)}`$ $``$ $`\{\begin{array}{cc}1\mathrm{GeV}\hfill & \mathrm{for}T_{RH}<T_{},\hfill \\ & \\ \mathrm{Min}[1\mathrm{GeV},f_{(2)}^{2n4}M_G^{n2}M^{(n3)}]\hfill & \mathrm{for}T_{RH}>T_{},M>f_{(1)}^{\frac{2}{n3}}f_{(2)}^{\frac{2(n4)}{n3}}M_G,\hfill \end{array}`$ (28)
$`M_F^{(u)}`$ $``$ $`\{\begin{array}{cc}\mathrm{Min}[\sqrt{m_{3/2}M_G},\left(m_\varphi ^{n1}M^{n3}\right)^{\frac{1}{2n4}}]\hfill & \mathrm{for}T_{RH}<T_{},\hfill \\ & \\ \mathrm{Min}[\sqrt{m_{3/2}M_G},\left(m_\varphi ^{n1}M^{n3}\right)^{\frac{1}{2n4}},f_{(2)}^{n1}M_G^{\frac{n1}{2}}M^{\frac{(n3)}{2}}]\hfill & \mathrm{for}T_{RH}>T_{},M>f_{(1)}^{\frac{2}{n3}}f_{(2)}^{\frac{2(n4)}{n3}}M_G,\hfill \end{array}`$ (32)
$`T_{RH}^{(l)}`$ $``$ $`f_{(2)}^2M_G^{\frac{1}{2}}\times \mathrm{Max}[\left(m_{3/2}^nM^{2(n3)}\right)^{\frac{1}{2(n2)}},\left(M_F^{2n}M^{n3}\right)^{\frac{1}{2(n1)}}],`$ (33)
$`T_{RH}^{(u)}`$ $``$ $`\{\begin{array}{cc}T_{inf}\hfill & \mathrm{for}M>f_{(1)}^{\frac{2}{n3}}f_{(2)}^{\frac{2(n4)}{n3}}M_G\hfill \\ & \\ \mathrm{Min}[f_{(1)}^{\frac{n}{n3}}f_{(2)}^{\frac{n+6}{n3}}M_G^{\frac{1}{2}}M^{\frac{3}{2}},T_{inf}]\hfill & \mathrm{for}M<f_{(1)}^{\frac{2}{n3}}f_{(2)}^{\frac{2(n4)}{n3}}M_G,\hfill \end{array}`$ (37)
where we do not take account of the gravitino problem, and the bounds for $`M_F`$ and $`m_{3/2}`$ include the constraints due to the scheme of the GMSB. Notice that $`T_{RH}^{(l)}`$ comes from $`\varphi _c>\mathrm{Max}[\varphi _\alpha ,\varphi _\gamma ]`$ and $`T_{RH}^{(u)}`$ is obtained from the condition $`f_{(1)}\varphi _c>T_c`$. Also, we must impose the following condition in order to take the gravitino problem Gouvea  into consideration.

$$T_{RH}<\gamma \times \{\begin{array}{cc}100\mathrm{G}\mathrm{e}\mathrm{V}1\mathrm{T}\mathrm{e}\mathrm{V}\hfill & \mathrm{for}2h^2\mathrm{keV}m_{3/2}100\mathrm{k}\mathrm{e}\mathrm{V}\hfill \\ & \\ 10\mathrm{T}\mathrm{e}\mathrm{V}\times h^2\left(\frac{m_{3/2}}{100\mathrm{k}\mathrm{e}\mathrm{V}}\right)\left(\frac{m_{G3}}{1\mathrm{T}\mathrm{e}\mathrm{V}}\right)^2\hfill & \mathrm{for}m_{3/2}>100\mathrm{k}\mathrm{e}\mathrm{V},\hfill \end{array}$$

where $`h`$ is the present Hubble parameter in units of 100km/sec/Mpc, $`m_{G3}`$ is the gaugino mass for $`SU(3)_C`$, and $`\gamma `$ reflects the entropy dilution due to the decay of Q-balls. If Q balls are subdominant component of the universe throughout, we have $`\gamma =1`$. On the other hand, if Q balls dominate the universe, it is given as

$$\gamma =\{\begin{array}{cc}\frac{T_{RH}}{T_d}\hfill & \mathrm{for}T_{eq}>T_{RH},\hfill \\ & \\ \frac{T_{eq}}{T_d}\hfill & \mathrm{for}T_{eq}<T_{RH},\hfill \end{array}$$

where $`T_{eq}`$ is the temperature at which the energy of Q balls becomes comparable to the total energy of the universe, and $`T_d`$ is the decay temperature of Q balls given as

$$T_d=\left(\frac{\pi ^2g_{}}{90}\right)^{\frac{1}{4}}\sqrt{\frac{M_FM_G}{48\pi Q^{\frac{5}{4}}}},$$
(38)
where $`g_{}`$ counts the effective degrees of freedom of the radiation.

So far we have outlined the dynamics of the amplitude of the AD field. Now we would like to focus on the motion of its phase, $`\theta \mathrm{arg}[\mathrm{\Phi }]`$, which essentially determines the baryon asymmetry. The question we have to ask here is when the baryon number is fixed. For this purpose let us begin our analysis by following the evolution of the baryon number density. For a moment, we do not consider the effect of the Q-ball formation. The equations of motion for $`\varphi `$ and $`\theta `$ are

$`\ddot{\varphi }+3H\dot{\varphi }+{\displaystyle \frac{V}{\varphi }}\dot{\theta }\varphi ^2`$ $`=`$ $`0,`$ (39)
$`\ddot{\theta }+\left(3H+2{\displaystyle \frac{\dot{\varphi }}{\varphi }}\right)\dot{\theta }+\varphi ^2{\displaystyle \frac{V}{\theta }}`$ $`=`$ $`0.`$ (40)
The baryon number density $`n_B`$ is given as

$`n_B`$ $`=`$ $`iq\left(\dot{\mathrm{\Phi }}^{}\mathrm{\Phi }\mathrm{\Phi }^{}\dot{\mathrm{\Phi }}\right),`$ (41)
$`=`$ $`q\dot{\theta }\varphi ^2,`$
where we defined that the baryon charge of the AD field is $`q`$. Its evolution is governed by the following equation,

$$\dot{n_B}+3Hn_B+q\frac{V}{\theta }=0.$$
(42)
The first and second terms mean that the number density is diluted by the expansion of the universe, and it is the third term that represents the generation of the baryon number. Note that the only A-term, $`V_A(+V_{AH})`$, can contribute to the third term.

In the first place we review how the baryon number is fixed in the ordinary AD mechanism. In this case, the AD field starts the oscillation when the Hubble parameter becomes equal to the mass scale of the radial direction, $`M_\varphi `$. Thereafter the amplitude $`\varphi `$ is reduced by the Hubble frictional term. Accordingly the magnitude of the A-term decreases rapidly, since it is a positive power of $`\varphi `$. Then it can be shown that the A-term cannot generate the baryon number efficiently after the AD field starts the oscillation as follows. While the AD field tracks the slowly decreasing minimum given as Eq. (15), the motion of $`\theta `$ is strongly damped by the frictional term because $`\dot{\varphi }/\varphi H`$ and $`HM_\varphi M_\theta `$, where $`M_\theta `$ is the mass scale of $`\theta `$. The second equality is satisfied if the Hubble-induced A-term exists. Hence the characteristic time scale for $`\theta `$ is $`(M_\theta ^2/H)^1`$, and the baryon number density is expressed as

$$n_B=q\dot{\theta }\varphi ^2q\frac{M_\theta ^2}{H}\varphi ^2.$$
(43)
Then the evolution equation for $`n_B`$ is given as

$`\dot{n_B}+3Hn_B`$ $`=`$ $`q{\displaystyle \frac{V}{\theta }},`$ (44)
$``$ $`qM_\theta ^2\varphi ^2,`$
$``$ $`Hn_B,`$
where we used Eq. (43) in the last equality. From Eq.(44), we can see that the baryon number continues to be critically generated while the AD field is trapped by the negative Hubble mass term. That is to say, the decrease of the baryon number density due to the expansion of the universe is compensated by the baryon number violating term to some extent, as long as the AD field tracks the instantaneous minimum given as Eq. (15). However, once it starts the oscillation, the magnitude of the A-term is significantly reduced, and also the time scale for $`\theta `$ becomes much smaller : $`\dot{\theta }^1M_\varphi ^1H^1`$. Thus the generation of the baryon number stops soon after the AD field starts oscillating in the ordinary AD mechanism.

The usual picture outlined above does not apply to the case when the AD field is trapped by the negative thermal log potential. First we assume the existence of the Hubble-induced A-term, and see how it spoils the baryogenesis. The relations among the mass scales, $`M_\varphi `$, $`M_\theta `$ and $`H`$, are then $`M_\varphi M_\theta H`$ before the trapping due to the thermal effect. Once the AD field is trapped by the negative thermal log potential, its amplitude, $`\varphi (f_{(2)}^2T^2M^{n3})^{1/(n1)}`$, decreases relatively slowly. While both $`M_\varphi `$ and $`M_\theta `$ are decreasing functions of slowly varying $`\varphi `$, the Hubble parameter diminishes faster. Thus, $`M_\varphi M_\theta H`$ holds until the trap ends. Since the mass scale of $`\theta `$ is larger than the Hubble parameter, $`\theta `$ starts oscillating around its minimum which is determined by the superposition of $`V_A`$ and $`V_{AH}`$, as soon as the AD field is trapped. During the trap, the amplitude of the oscillation of $`\theta `$ decreases continuously 222Though the Hubble-induced A-term, $`V_{AH}`$, might become smaller than the A-term, $`V_A`$, it does not increase the amplitude of $`\theta `$., so that $`\theta `$ stays very close to the bottom of the valley of the A-term when $`\varphi `$ starts an oscillation. Therefore the AD field rotates the orbit with very large oblateness, as shown schematically in Fig. 2. It should be noted that the amplitude of the oscillation of $`\varphi `$ remains almost constant due to the small Hubble parameter, which has a decisive influence on the evolution of the baryon number density. In order to follow the evolution of $`n_B`$, we estimate the third term in Eq. (42), which is now given as

$`q{\displaystyle \frac{V}{\theta }}`$ $``$ $`qM_\theta ^2\varphi ^2=\dot{\theta }^1M_\theta ^2n_B,`$ (45)
$``$ $`M_\theta n_B,`$
where the last equality holds at the end point of each oscillation. Therefore the baryon number density continues to be generated around the end point, since $`M_\theta H`$. Put simply, the valleys of the A-term become effectively deeper due to the enduring trap, which causes the intermittent generation of the baryon number.

Consider now the implication of the Q-ball formation. The AD field experiences instabilities when it oscillates in the $`V_{gauge}`$ dominated potential, and deforms into the gauge-mediation type Q balls. Since the growth rate of the instability is of the order of $`M_\varphi M_F^2/\varphi _{osc}`$ Kasuya3 , Q balls are formed soon after the AD field starts a radial oscillation 333The existence of the A-term does not prevent the formation of Q balls, because it does not alter the flatness of the potential at all, namely $`\mathrm{Max}[V_A,V_{AH}]<V_{gauge}`$. In other words, if only a small part of the energy for the radial component is converted to that for the phase direction, the AD field can easily get over the hill of the A-term. Also we confirmed it with the numerical calculation.. The most important part of this argument is that the baryon number is not fixed during the process toward the Q-ball formation. When the fluctuation becomes comparable to the homogeneous mode, the system goes into the nonlinear stage, and the AD field is kicked back and forth at random in the phase direction by the A-term. It is not until Q balls are formed that the baryon number is fixed. The reason is as follows. It is known that Q balls absorb almost all the baryon number. Inside Q balls, the AD field rotates a circular orbit, so the average of the third term in Eq. (42) over a round of the orbit vanishes:

$$_0^{2\pi }𝑑\theta \frac{V}{\theta }=0,$$
(46)
which means the baryon number is fixed. Since the baryon number evolves indiscriminately until the Q-ball formation is completed, the resultant baryon asymmetry is affected by the highly nonlinear physics of the Q-ball formation. In rough estimate, the created baryon number before the system goes into the nonlinear stage is

$$n_B^{(i)}\frac{c_1}{M_\varphi }M_\theta ^2\varphi _{osc}^2\mathrm{sin}n\theta _{osc},$$
(47)
where $`\theta _{osc}(1)`$ is the phase amplitude at the start of the radial oscillation, and $`c_1M_\varphi ^1`$ with $`c_1\mathrm{log}(\varphi /\delta \varphi )=O(10)`$ is the time for the fluctuation to become nonlinear. When the fluctuation of the phase becomes large, i.e., $`\delta \theta /\theta O(1)`$, the AD field moves beyond the potential hill in the phase direction and the baryon number production restarts with random phase, $`\theta _i(𝐱)O(1)`$. Here $`\theta _i`$ depends on the position $`𝐱`$. Hence this process induces the fluctuation of the baryon number density, and it takes place over and over again until Q-ball formation is completed. If the process lasted eternally, however large the baryon asymmetry is, it would be completely erased due to the baryon number violating process in equilibrium. Note that the baryon number is odd under the CPT transformation. In fact, this process is effective for a finite time, so that the only finite baryon asymmetry is annihilated. We can estimate the magnitude of the baryon asymmetry, which is erased due to the Q-ball formation processes, by calculating the typical value of the fluctuation of the baryon number density :

$`\delta n_B(𝐱)`$ $``$ $`{\displaystyle \frac{1}{M_\varphi }}M_\theta ^2\varphi _{osc}^2{\displaystyle \underset{i}{\overset{c_2}{}}}\mathrm{sin}n\theta _i(𝐱),`$ (48)
$``$ $`\pm {\displaystyle \frac{\sqrt{c_2}}{M_\varphi }}M_\theta ^2\varphi _{osc}^2,`$
where $`c_2O(10)`$ is the number of the restarts, and the contribution of the random phase $`\theta _i`$ is evaluated as random walk with $`c_2`$ steps in the second line. Since $`\theta _{osc}1`$ in the present situation, we have $`n_B^{(i)}\delta n_B`$, which means that $`n_B`$ is dominated by the fluctuation due to the Q-ball formation process. Moreover, the dynamics is chaotic in the sense that a tiny change results in a quite different baryon number. Thus, the baryon asymmetry generated during the homogeneous evolution is swallowed up in the large fluctuation, leading to the baryon symmetric universe.

We performed the numerical simulation in order to confirm whether the baryon asymmetry is really determined through the Q-ball formation. For the numerical calculation, we take the variables dimensionless as follows.

$`\phi `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }}{M_S}},`$
$`\tau `$ $`=`$ $`{\displaystyle \frac{M_F^2}{M_S}}t,`$
$`\chi _i`$ $`=`$ $`{\displaystyle \frac{M_F^2}{M_S}}x_i,`$ (49)
The initial conditions of the homogeneous parts are taken as

$`\mathrm{Re}\left[\phi (0)\right]`$ $`=`$ $`1.0\times 10^2,`$
$`\mathrm{Im}\left[\phi (0)\right]`$ $`=`$ $`1.0\times 10^3,`$
$`\phi ^{}(0)`$ $`=`$ $`0.0,`$ (50)
where the prime denotes the derivative with respect to $`\tau `$, and we redefined the phase of the AD field so that the real direction coincides with one of the valleys of the A-term with the mass scale $`M_\theta /M_\varphi =10^{1.8}`$. Taking account of fluctuations which originate from the quantum fluctuations, we give $`\phi `$ fluctuations with amplitudes $`10^7`$ times smaller than the homogeneous mode. We have confirmed that the smaller fluctuations just delay the formation of Q balls. For simplicity we neglected the cosmic expansion. The results are shown in Fig. 3. It shows that the baryon number is generated at the end point of each oscillation for some time $`c_1M_\varphi ^1`$, but it soon goes into the chaotic phase until the Q-ball formation is completed. If we change the initial seed for the fluctuation, it is probable that the resultant baryon asymmetry inside the box of the lattices becomes positive. The reason why the resultant baryon number takes nonzero value is that the size of the box under consideration is finite. If the same calculation were done with the infinitely large box, the resultant baryon number would be exactly zero.

Lastly we comment on the case without the Hubble-induced A-term. Now $`M_\varphi >H>M_\theta `$ at the onset of the oscillation can be realized since the mass scale of $`\theta `$ is smaller than that in the previous case. Then the initial phase amplitude $`\theta _{osc}`$ is comparable to the random phase $`\theta _i`$, hence the net baryon number can be nonzero. In the next section we estimate the resultant baryon asymmetry in this case, taking account of the Q-ball formation.

## IV Baryogenesis

As we have shown in the previous section, the Hubble-induced A-term spoils the baryogenesis through the Q-ball formation. On the other hand, if it were not for the Hubble-induced A-term, we can avoid this problem since $`M_\varphi >H>M_\theta `$ at the beginning of the oscillation is possible. Then the baryon number just before the Q-ball formation is completed is roughly estimated as

$$n_B^{(f)}\frac{c_1}{M_\varphi }M_\theta ^2\varphi _{osc}^2>\frac{\sqrt{c_2}}{M_\varphi }M_\theta ^2\varphi _{osc}^2,$$
(51)
where $`c_1c_2O(10)`$, and we used $`\theta _{osc}O(1)`$ in the present case. Thus we expect that the finite baryon asymmetry remains even if the AD field gets randomly kicked in the phase direction during the nonlinear stage. In other words, the resultant baryon number is almost determined due to the dynamics before the AD field goes into the nonlinear phase, so subsequent Q-ball formation does not affect the baryon asymmetry.

To satisfy the requirement that the net baryon number is generated, the following condition must be satisfied:

$$M_\theta <H_{osc},$$
(52)
where

$`M_\theta `$ $``$ $`\sqrt{{\displaystyle \frac{m_{3/2}\varphi _{osc}^{n2}}{M^{n3}}}},`$ (53)
$`H_{osc}`$ $``$ $`\{\begin{array}{cc}{\displaystyle \frac{T_{osc}^4}{T_{RH}^2M_G}}\hfill & \mathrm{for}T_{osc}>T_{RH},\hfill \\ & \\ {\displaystyle \frac{T_{osc}^2}{M_G}}\hfill & \mathrm{for}T_{osc}<T_{RH}.\hfill \end{array}`$ (57)
Then the baryon number is generated as

$$n_B(t_{osc})\frac{1}{M_\varphi }\frac{m_{3/2}\varphi _{osc}^n}{M^{n3}},$$
(58)
where it should be noted that $`M_\varphi =M_F^2/\varphi _{osc}`$ is used instead of the Hubble parameter, because the kick due to the A-term lasts for $`\delta tM_\varphi ^1`$, and we also set the numerical factor to be unity. The almost all baryon number is absorbed into the subsequently formed Q balls with the typical charge given as Eq. (12444 The value of $`\beta `$ in Eq.(12) might be different in our case, since the Q-ball formation is delayed compared to the ordinary scenario. However we make use of this equation, since it is still valid up to the numerical factor.. Hereafter we consider the two cases, depending on whether the Q balls dominate the universe or not. If the Q balls dominates the universe, they must be unstable and decay before the epoch of big bang nucleosynthesis (BBN). Thus, the mass per unit charge and the decay temperature are constrained as

$`M_FQ^{\frac{1}{4}}`$ $`>`$ $`1\mathrm{G}\mathrm{e}\mathrm{V},`$ (59)
$`T_d=\left({\displaystyle \frac{\pi ^2g_{}}{90}}\right)^{\frac{1}{4}}\sqrt{{\displaystyle \frac{M_FM_G}{48\pi Q^{\frac{5}{4}}}}}`$ $`>`$ $`10\mathrm{M}\mathrm{e}\mathrm{V}.`$ (60)
Since the universe is reheated again, the gravitino problem is alleviated in this case. Accordingly the constraints Eqs.(III) and (III) with $`\gamma `$ given as Eq. (III) are adopted. The baryon-to-entropy ratio is given as

$`{\displaystyle \frac{n_B}{s}}`$ $`=`$ $`{\displaystyle \frac{3}{4}}{\displaystyle \frac{n_B(t_{osc})T_d}{M_F^4}},`$ (64)
$``$ $`\{\begin{array}{cc}30\times \left({\displaystyle \frac{g_{}}{10}}\right)^{\frac{1}{4}}m_{3/2}^{n+\frac{5}{2}}M_F^{2(n3)}M^{(n3)}M_G^{\frac{1}{2}}\hfill & \mathrm{for}m_{3/2}>\left(M_F^{2n4}M^{(n3)}\right)^{\frac{1}{n1}},\hfill \\ & \\ 30\times \left({\displaystyle \frac{g_{}}{10}}\right)^{\frac{1}{4}}m_{3/2}M_F^{\frac{n}{n1}}M^{\frac{(n3)}{2(n1)}}M_G^{\frac{1}{2}}\hfill & \mathrm{for}m_{3/2}<\left(M_F^{2n4}M^{(n3)}\right)^{\frac{1}{n1}},\hfill \end{array}`$
The allowed regions for $`m_{3/2}`$, $`M_F`$ and $`T_{RH}`$ satisfying the above constraints are shown in Fig.4 and 5, where we set $`f_{(1)}=f_{(2)}=0.3`$, $`M=M_G`$ and $`n=6`$. It shows that the baryon-to-entropy ratio can be large enough to explain the present asymmetry. For instance, it takes the value as

$$\frac{n_B}{s}=6\times 10^7\left(\frac{g_{}}{10}\right)^{\frac{1}{4}}\left(\frac{m_{3/2}}{100\mathrm{M}\mathrm{e}\mathrm{V}}\right)\left(\frac{M_F}{5\times 10^8\mathrm{GeV}}\right)^{\frac{6}{5}}\left(\frac{M}{M_G}\right)^{\frac{3}{10}},$$
(65)
where $`m_{3/2}`$ and $`M_F`$ is assumed to satisfy all the constraints stated above, so their values must be included in the allowed regions shown in Fig.4 and 5. Though the above value is a little larger compared to the expected one $`10^{10}`$, it might be suppressed by the unknown CP phase $`\theta _{osc}`$which we assume to be of the order of unity here. Although $`M>M_G`$ is possible since the definition of $`M`$ in Eq. (3) includes a possibly small coupling constant, it does not work to take larger $`M`$. This is because the requirement that Q balls decay before the BBN sets the lower bound for $`M_F`$ which is a positive power of $`M`$, so the allowed regions shown in Fig.4 and 5 would disappear for $`MM_G`$.

Next we consider the case that the energy of Q balls is subdominant before the BBN. If Q balls are unstable, they must decay before the BBN, otherwise their decay product might destroy the light elements and spoil the success of the BBN. If Q balls are stable, the resultant baryon asymmetry is suppressed by $`\mathrm{\Delta }Q/Q`$, where $`\mathrm{\Delta }Q`$ is the evaporated charge. Also the abundance of remnant Q balls are constrained in order not to overclose the universe. Namely,

$$\mathrm{\Omega }_Q=\frac{M_FQ^{1/4}}{m_n}\frac{\mathrm{\Omega }_B}{\epsilon }\frac{Q}{\mathrm{\Delta }Q}<1,$$
(66)
where $`m_n`$ is the nucleon mass, and $`\mathrm{\Omega }_{Q(B)}`$ is the ratio of the energy density of Q balls (or baryon) to the critical density. The baryon-to-entropy ratio is expressed as

$$\frac{n_B}{s}\left(\frac{45}{2\pi ^2g_{}}\right)\frac{n_B(t_{osc})}{T_{osc}^3}\times \kappa _1\times \kappa _2,$$
(67)
where

$`\kappa _1`$ $``$ $`\{\begin{array}{cc}1\hfill & \mathrm{for}T_{osc}<T_{RH},\hfill \\ & \\ T_{RH}^5/T_{osc}^5\hfill & \mathrm{for}T_{osc}>T_{RH},\hfill \end{array}`$ (71)
$`\kappa _2`$ $``$ $`\{\begin{array}{cc}1\hfill & \mathrm{for}\mathrm{unstable}\mathrm{Q}\mathrm{balls},\hfill \\ & \\ \mathrm{\Delta }Q/Q\hfill & \mathrm{for}\mathrm{stable}\mathrm{Q}\mathrm{balls}.\hfill \end{array}`$ (75)
Since the reheating temperature is now constrained by the gravitino problem, the constraints Eqs.(III) and (III) with $`\gamma =1`$ are adopted. The allowed regions for $`m_{3/2}`$, $`M_F`$ and $`T_{RH}`$ satisfying the above constraints are shown in Fig.6 and 7, where we set $`f_{(1)}=f_{(2)}=0.3`$, $`M=M_G`$ and $`n=6`$. It shows that the baryon-to-entropy ratio can be as large as $`10^{12}`$, which is a little smaller than the favored value $`10^{10}`$. However, it can be cured by taking smaller $`M`$, the GUT scale for instance. Fig.8 and 9 represent the allowed regions for $`m_{3/2}`$, $`M_F`$ and $`T_{RH}`$ in the case of $`f_{(1)}=f_{(2)}=0.3`$, $`M=10^{16}`$ GeV and $`n=6`$. It shows that the baryon-to-entropy ratio can be now as large as $`10^{10}`$ for unstable Q balls, although it is too small to explain the present asymmetry for stable Q balls. For instance, it is given as

$$\frac{n_B}{s}3\times 10^{11}\left(\frac{g_{}}{200}\right)^1\left(\frac{m_{3/2}}{100\mathrm{M}\mathrm{e}\mathrm{V}}\right)\left(\frac{M_F}{2\times 10^7\mathrm{GeV}}\right)^{\frac{36}{5}}\left(\frac{T_{RH}}{10^6\mathrm{GeV}}\right)^5\left(\frac{M}{10^{16}\mathrm{GeV}}\right)^{\frac{6}{5}},$$
(76)
for unstable Q balls.

## V Conclusion

We have studied whether the AD baryogenesis works in the gauge-mediated SUSY breaking models, if the AD field is trapped by the negative thermal log potential which is known to be ubiquitous in the MSSM. The most striking feature is that the Hubble parameter becomes much smaller than the mass scale of the radial component of the AD field, during the trap. Then, the amplitude of the AD field decreases so slowly that the baryon number is not fixed even after the onset of radial oscillation. A further important point is that the AD field experiences spatial instabilities and deforms into Q balls. Therefore the baryon asymmetry can be washed out due to the random kick in the phase direction during the process of Q-ball formation. Indeed, it has been shown that the net baryon number vanishes if the Hubble parameter is smaller than the mass scale of phase direction, since the randomly generated baryon number surpasses the baryon number just before the nonlinear stage begins. This is always the case if the Hubble-induced A-term exists. Also we confirmed the chaotic behavior of the baryon number with use of numerical calculation. However, the Hubble-induced A-term does not necessarily exist. In fact, we have found that it is possible to generate the right amount of the baryon asymmetry in the absence of the Hubble-induced A-term, if the Q balls are unstable. Thus, it should be emphasized that the AD baryogenesis is still a viable candidate, even if both the thermal effect and Q-ball formation are fully taken into account.

### ACKNOWLEDGMENTS

F.T. thanks Masahide Yamaguchi for useful discussion and comments.