text:
# Bosonic model with 𝑍₃ fractionalization

## I Introduction

A flurry of recent theoretical activity has produced specific model system realizations of fractionalized phases in two dimensions.RSSpN ; MoeSon ; Iof ; BalMPAFGir ; frcmdl ; bosfrc  Essentially all of the fractionalized states constructed so far are $`Z_2`$ states.U1  On a formal level, these realizations employ the following route to $`Z_2`$ fractionalization: Strong local correlations lead to a $`U(1)`$ gauge theory as a low-energy description; this gauge theory is then driven into a deconfined state by a condensation of objects carrying gauge charge $`2`$ that also appear in the low-energy description. This formal structure has been brought out very directly in Refs. SV, ; frcmdl, ; bosfrc, . On a more physical level, the fractionalized insulator is produced departing from a superconducting state by a condensation of double vortices.NLII ; z2long  The main body of work concentrated on the $`Z_2`$ states since these are expected to be the simplest to realize. However, it is clear that more complicated fractionalized states are also possible. For example, it is conceivable that in some system the superconducting state is quantum-disordered by a condensation of triple vortices; the resulting insulator is then a $`Z_3`$ fractionalized state.

In this paper, we indicate how a $`Z_3`$ fractionalized state can be engineered in a relatively simple bosonic model with unfrustrated nearest-neighbor hopping and short-range two-body repulsive interaction. Much of the construction parallels closely the $`Z_2`$ examples of Refs. frcmdl, ; bosfrc, ; Iof, : The low-energy Hilbert space is selected—by stipulating particular charge interactions—in a manner that naturally admits splitting boson charge into three pieces; this Hilbert space is protected by a large charge gap. The effective description of the fractionalized state has gapped chargons carrying electrical charge $`1/3`$ and coupled to some special $`Z_3`$ gauge theory which we analyze in detail. Our main message here is that one does not need very contrived systems to obtain more complicated fractionalization patterns.

## II $`Z_3`$ via charge frustration

The model is defined on the lattice shown in Fig. 1, which we can think of as a hexagonal lattice with additional sites placed at the hexagon centers. We have $`\psi _r^{}=e^{i\varphi _r}`$ bosons residing on the hexagonal lattice (these sites are always labelled lower-case $`r`$), and $`b_R^{}=e^{i\theta _R}`$ bosons residing at the hexagon centers (upper-case $`R`$). Bosons can hop between the neighboring sites as indicated by the links on the figure; the hopping amplitudes are $`w_1`$ for $`r,R`$ links and $`w_2`$ for $`r,r^{}`$ links. We also stipulate strong repulsive interactions that favor charge neutrality of the hexagons, in addition to the on-site repulsion that favors charge neutrality of the individual sites. The complete quantum rotor Hamiltonian is

$`H`$ $`=`$ $`w_1{\displaystyle \underset{R,rR}{}}(b_R^{}\psi _r+\mathrm{h}.\mathrm{c}.)w_2{\displaystyle \underset{rr^{}}{}}(\psi _r^{}\psi _r^{}+\mathrm{h}.\mathrm{c}.)`$ (1)
$`+u_b{\displaystyle \underset{R}{}}(n_R^b)^2+u_\psi {\displaystyle \underset{r}{}}(n_r^\psi )^2+U{\displaystyle \underset{R}{}}N_R^2.`$
Here, $`\{n_R^b,\theta _R\}`$ are conjugate number-phase variables \[e.g., in the phase representation $`n_R^bi/\theta _R`$, $`b_R^{}\psi _r+\mathrm{h}.\mathrm{c}.2\mathrm{cos}(\theta _R\varphi _r)`$\], and similarly for $`\{n_r^\psi ,\varphi _r\}`$. The number-phase variables are particularly appropriate if we think of the model as describing an array of Josephson-coupled superconducting islands. Both $`b_R`$ and $`\psi _r`$ bosons carry electrical charge $`q_b`$.

In the above Hamiltonian, $`N_R`$ is the boson number associated with each hexagon:

$$N_R=3n_R^b+\underset{rR}{}n_r^\psi .$$
(2)
Thus, the total boson number in the systemboundaries  is

$$N_{\mathrm{tot}}=\frac{1}{3}\underset{R}{}N_R.$$
(3)
General analysis of the possible phases proceeds as in Ref. bosfrc, . Here, we focus on the insulating states that are obtained for large $`Uw_1,w_2,u_b,u_\psi `$. If the $`w_1,w_2,u_b,u_\psi `$ terms are all zero, there is a degenerate manifold of ground states specified by the requirement $`N_R=0`$ for each $`R`$.numbers  The ground state sector is separated by a large charge gap $`U`$ from the nearest sectors. Including the $`w_1,w_2,u_b,u_\psi `$ terms lifts this degeneracy in each sector, which is best described by deriving the corresponding effective Hamiltonian in each sector.

The effective Hamiltonian in the ground state sector ($`N_R=0`$) is, to third order in the perturbing terms,

$`H_{\mathrm{eff}}^{(0)}=H_{u_b,u_\psi }J_\mathrm{c}{\displaystyle \underset{r}{}}[(\psi _r^{})^3b_{R_1}b_{R_2}b_{R_3}+\mathrm{h}.\mathrm{c}.]`$
$`K_{\mathrm{ring}}{\displaystyle \underset{\mathrm{}}{}}(\psi _1^{}\psi _2\psi _3^{}\psi _4\psi _5^{}\psi _6+\mathrm{h}.\mathrm{c}.).`$ (4)
Here, $`H_{u_b,u_\psi }`$ stands for the on-site repulsion terms as in Eq. (1); $`R_1,R_2,R_3`$ label the three hexagon centers adjacent to $`r`$; $`J_\mathrm{c}=w_1^3/(6U^2)`$ and $`K_{\mathrm{ring}}=3w_2^3/U^2`$.

This is our main step in obtaining $`Z_3`$ fractionalization. The above Hamiltonian looks similar to a compact $`U(1)`$ gauge theory coupled to a charge $`3`$ scalar field. Thus, if we think of the $`\psi _r`$ as some gauge fields, then it is very suggestive to think of the $`b_R`$ as carrying gauge charge $`3`$ \[see also Eq. (2)\]. Standard Fradkin-Shenker analysisFraShe  then suggests that by condensing the $`b_R`$ field, which can be arranged by making $`J_\mathrm{c},K_{\mathrm{ring}}`$ large, we can deconfine objects carrying gauge charge $`1`$. By virtue of Eq. (3), such objects carry fractional electrical charge $`q_b/3`$, and we obtain a $`Z_3`$ fractionalized insulator.

The reader who finds the above statements believable may now declare victory in achieving $`Z_3`$ fractionalization. However, if we want to describe the deconfined phase(s) in any detail, we need to study the above Hamiltonian directly since it is not related in any simple manner to the conventional gauge theory. This is our focus in the remainder of the paper. Of course, we will confirm the deconfinement, but we will also find that the deconfined state that obtains for large $`K_{\mathrm{ring}}`$ on all hexagons is in fact a $`Z_3\times Z_3`$ state (see below for details).

Proceeding with this analysis, consider the regime of large $`J_\mathrm{c}\mathrm{}`$ and small $`u_b0`$. It is convenient to perform the following change of variables. Define the operators $`b_{cR}^{}=e^{i\theta _{cR}}`$ and $`\stackrel{~}{\psi }_r^{}=e^{i\stackrel{~}{\varphi }_r}`$:

$$b_{cR}^{}=s_Re^{i\theta _R/3},\stackrel{~}{\psi }_r^{}=\psi _r^{}b_{cR_1}b_{cR_2}b_{cR_3}.$$
(5)
Here, $`s_R=1,e^{i2\pi /3}`$, or $`e^{i4\pi /3}`$, so that $`\theta _{cR}[0,2\pi )`$. One can readily verify that $`N_R`$ is conjugate to $`\theta _{cR}`$, while $`n_r^\psi `$ is conjugate to $`\stackrel{~}{\varphi }_r`$. $`b_{cR}^{}`$ carries electrical charge $`q_b/3`$ and can be thought of as a chargon field, while $`\stackrel{~}{\psi }_r^{}`$ is charge neutral. These new variables are indeed natural in the description of the deconfined phase, but to recover the physical Hilbert space, we need to impose the constraint

$`\mathrm{exp}\left[i{\displaystyle \frac{2\pi }{3}}(N_R{\displaystyle \underset{rR}{}}n_r^\psi )\right]=1.`$ (6)
In the new variables, the Hamiltonian becomes

$`H_{\mathrm{eff}}^{(0)}=u_\psi {\displaystyle \underset{r}{}}(n_r^\psi )^2J_\mathrm{c}{\displaystyle \underset{r}{}}[(\stackrel{~}{\psi }_r^{})^3+\mathrm{h}.\mathrm{c}.]`$
$`K_{\mathrm{ring}}{\displaystyle \underset{\mathrm{}}{}}(\stackrel{~}{\psi }_1^{}\stackrel{~}{\psi }_2\stackrel{~}{\psi }_3^{}\stackrel{~}{\psi }_4\stackrel{~}{\psi }_5^{}\stackrel{~}{\psi }_6+\mathrm{h}.\mathrm{c}.).`$ (7)
Note that we are left with gauge fields $`\stackrel{~}{\psi }_r`$ only, since we are working in the “uncharged” ground state sector

The term $`J_\mathrm{c}`$ acts as a $`Z_3`$ anisotropy on the $`\stackrel{~}{\varphi }_r`$ field. In the limit $`J_\mathrm{c}\mathrm{}`$, $`e^{i\stackrel{~}{\varphi }_r}`$ becomes a $`Z_3`$ field: $`\stackrel{~}{\varphi }_r=0,2\pi /3`$, or $`4\pi /3`$. The operator $`P_r^+e^{i(2\pi /3)n_r^\psi }`$ shifts the states of the quantum $`Z_3`$ clock by $`+1`$,ZN  whereas the constraints specifying the uncharged sector $`N_R=0`$ become

$`{\displaystyle \underset{rR}{}}P_r^+=1.`$ (8)
The $`u_\psi (n_r^\psi )^2`$ term causes tunneling between the different states of the quantum clock, and this can be described by an effective “transverse field” $`h(P_r^++P_r)`$ acting on the clock.

The effective Hamiltonian now becomes $`Z_3`$ “ring-exchange” Hamiltonian on the hexagonal lattice of $`r`$ sites

$`K_{\mathrm{ring}}{\displaystyle \underset{\mathrm{}}{}}(\stackrel{~}{\psi }_1^{}\stackrel{~}{\psi }_2\stackrel{~}{\psi }_3^{}\stackrel{~}{\psi }_4\stackrel{~}{\psi }_5^{}\stackrel{~}{\psi }_6+\mathrm{h}.\mathrm{c}.)h{\displaystyle \underset{r}{}}(P_r^++P_r^{}).`$ (9)
The Hamiltonian together with the constraints Eq. (8) can be viewed as some special $`Z_3`$ gauge theory and is analyzed below and in further detail in Appendices A and B. We find that generically this theory can have two deconfined phases (in addition to the confined phase) with the phase diagram shown in Fig. 2. Here we only describe the deconfined phase that obtains when all ring exchange couplings are large, $`K_{\mathrm{ring}}h`$. As explained below, this phase is a $`Z_3\times Z_3`$ deconfined phase.

From here on, our focus is on the above $`Z_3`$ ring exchange Hamiltonian. We drop all tildes on $`Z_3`$ fields $`\psi _r`$ and superscripts on $`n_r`$ (which are now integers modulo $`3`$). Also, we consider a $`Z_N`$ generalizationZN  of the above Hamiltonian and carry out the analysis in the general case. This is done for clarity of notation.

For $`K_{\mathrm{ring}}h`$, a good caricature of the bulk ground state is given by the wavefunction

$`|GS={\displaystyle \underset{\{n_r\}}{\overset{}{}}}|\{n_r\},`$ (10)
where the primed sum is over all configurations $`\{n_r\}`$ that satisfy the constraints Eq. (8), i.e., $`_{rR}n_r=0`$.

Let us define $`Z_N`$ flux through a given hexagon $`R`$

$`\mathrm{\Phi }_R=\varphi _1\varphi _2+\varphi _3\varphi _4+\varphi _5\varphi _6,`$ (11)
with the sign convention as in Fig. 3. The ground state has zero flux through each hexagon. Excitations above this ground state are $`Z_N`$ vortices. For example, we can add a unit of flux through a given hexagon by applying a “string” operator as indicated in Fig. 3. The gap for a vortex carrying one unit of flux is $`2K_{\mathrm{ring}}[1\mathrm{cos}(2\pi /N)]`$.

Observe now (Fig. 3) that the hexagon centers $`R`$ form a triangular lattice, which consists of three sublattices A, B, and C. Observe also that the flux-adding string operator “steps” only through the same sublattice hexagons. We are thus led to the possibility of a topological distinction between vortices on the different sublattices, in addition to the usual distinction between two vortices carrying different flux. Indeed, one can see that the topologically distinct situations can be characterized by saying that we have two species of $`Z_N`$ vortices, say A and B vortices. Alternatively, if we want to preserve the symmetry among the three sublattices, we can say that there are three types of $`Z_N`$ vortices—A, B, and C vortices—but these are not independent and instead satisfy “fusion rules” such as

$`\left(\mathrm{\Phi }_A=+1\right)\times \left(\mathrm{\Phi }_B=+1\right)\left(\mathrm{\Phi }_C=1\right).`$ (12)
This means that a nearby pair of $`+1`$ A and B vortices is indistinguishable from a $`1`$ C vortex (note also that the “states” on the left and on the right can be connected by local $`h`$ terms in the Hamiltonian).

Consider now introducing $`Z_N`$ charges in the above gauge theory Eq. (9), e.g., consider placing a pair of opposite $`\pm 1`$ charges on two hexagons $`R_1`$ and $`R_2`$: $`_{rR}n_r=\delta _{RR_1}\delta _{RR_2}`$. This is appropriate when studying the charged sectors of the microscopic Hamiltonian Eq. (1) since the formal gauge structure represents the crucial coupling of chargons with the above $`Z_3`$ degrees of freedom. A charge can be added to a hexagon by applying a string operator as indicated in Fig. 3. From several perspectives, one can see that all such charges are deconfined in the $`K_{\mathrm{ring}}h`$ phase: Thus, in Appendix A we approach this “fully deconfined” phase starting from a “partially deconfined” phase, in which charges are deconfined on one sublattice only. Also, this fully deconfined phase corresponds to the fully disordered phase in the dual global $`Z_N`$ spin model of Appendix B.

Similarly to vortices, we need to distinguish the charges on different sublattices. Again, as far as the gauge structure is concerned, we have fusion rules such as

$`\left(N_A=+1\right)\times \left(N_B=+1\right)\left(N_C=1\right).`$ (13)
Statistical interactions between the different particles are readily identified by studying the commutation properties of the corresponding strings. These are summarized in Table 1 by specifying “gauge charges” of the different $`N_R=+1`$ excitations with respect to the A, B, and C fluxes. Thus, $`N_A=+1`$ excitation carries gauge charges $`Q_A=0,Q_B=1,Q_C=+1`$, i.e., it does not “see” A vortices, but when transported around a B or C vortex of unit strength, the wavefunction acquires an additional phase $`e^{i2\pi /N}`$ or $`e^{i2\pi /N}`$ correspondingly.

This completes the particle description of the fully deconfined phase. The minimal description would be to say that we have A hexagon charges that see B hexagon $`Z_N`$ vortices and B hexagon charges that see A hexagon vortices. This is essentially the claimed $`Z_N\times Z_N`$ structure. Thus, we expect $`N^2`$-fold ground state degeneracy if the system is put on a cylinder, as can be verified by constructing the corresponding ground states starting from the state Eq. (10) and threading vortices through the hole of the cylinder.

It should be emphasized here that the above discussion assumed that the three-sublattice structure is respected by the boundary conditions. While it is clear that the bulk properties do not depend on this, there is an additional quirk when we consider topological degeneracy in a geometry that does not respect the three-sublattice structure. This is legitimate when all ring exchange couplings are equal. Consider, e.g., a cylindrical geometry with the circumference along the horizontal direction of Figs. 1 and 3. When an A-type particle is transported around the periodic direction, it does not return to its initial position, but rather becomes a B\- or C-type particle. It takes three turns for the particle to return to the original position. From table 1, such $`N_R=+1`$ particle will not register any flux in this process. A detailed analysis shows that for $`N\mathrm{mod}\mathrm{\hspace{0.33em}3}0`$ the ground state of the system in this geometry is nondegenerate. On the other hand, for $`N\mathrm{mod}\mathrm{\hspace{0.33em}3}=0`$ the ground state is found to be three-fold degenerate, since in this case there is a composite object that returns to its initial state when transported once around the cylinder and that senses some flux through the hole of the cylinder in the process.

Returning to our microscopic bosonic model, the “particle description” of the $`Z_3\times Z_3`$ phase is as follows: We have two species of $`Z_3`$ vortices (with gap $`K_{\mathrm{ring}}`$) and we have charged particles (with charge gap $`U`$) that can be classified as carrying two distinct $`Z_3`$ gauge charges, in addition to their fractional electrical charge. Finally, note that the $`Z_3\times Z_3`$ state is associated with the additional symmetries in the hexagonal lattice ring exchange Hamiltonian but is protected by the same charge gap projection, since any move within the uncharged state sector is necessarily a combination of hexagon ring exchanges.

## III Conclusions

We showed that it is possible to produce more complicated fractionalization patterns such as $`Z_3`$ fractionalization in relatively simple bosonic models. While the resulting fractionalized state turned out to be even more complicated than initially intended, the microscopic model was not too contrived. It is hoped that this work will encourage further searches for other exotic states. For example, can a non-Abelian fractionalized state be produced in a condensed matter system with a global symmetry only, short-range interactions, and in zero magnetic field?

###### Acknowledgements.

The author particularly thanks T. Senthil for many stimulating discussions and persistent encouragement that this toy project be written up, and also for useful comments on the manuscript. This work was supported by the MRSEC program of the National Science Foundation under grants DMR-9808941 and DMR-0213282.

## Appendix A Partially deconfined ($`Z_N`$) phase

To better appreciate the character of the deconfinement in the special gauge theory Eq. (9), we allow different ring exchange couplings for different hexagons and consider particular parameter space with two such couplings: $`K_{\mathrm{ring}}=K_A`$ for the A hexagons and $`K_{\mathrm{ring}}=K_B=K_CK`$ for the B and C hexagons. This is indicated schematically in Fig. 3 where the A hexagons are shaded. Note that by allowing the two couplings we implicitly assume that the boundary conditions on the lattice respect the three-sublattice structure; this is done throughout.

We argue below that the ring exchange Hamiltonian has the phase diagram shown in Fig. 2 with three phases: For $`hK,K_A`$ the system is in a fully confined phase. For $`Kh\sqrt{K_AK}`$ the system is in a partially deconfined ($`Z_N`$) phase. In this phase, the charges on the A hexagons are deconfined, while the charges on the B and C hexagons are confined. Finally, for $`K,K_Ah`$ the system is in a fully deconfined ($`Z_N\times Z_N`$) phase with all charges deconfined. The phase diagram of Fig. 2 is also supported by the analysis of the dual global $`Z_N`$ spin model summarized in Appendix B.

In what follows, we give a detailed description of the partially deconfined phase. As a representative of this phase, consider the Hamiltonian with $`K_A=0`$, i.e., with ring exchanges around the B and C hexagons only (see Fig.3). In this case, there are additional conserved quantities:

$`\widehat{}_{AA^{}}n_r+n_r^{}=\mathrm{const}\text{(model with }K_A=0\text{)}`$ (14)
for each hexagonal lattice link $`rr^{}`$ between two A hexagons $`A`$ and $`A^{}`$ (see Fig. 3). This facilitates the analysis, since we can consider separately each subsector specified by the corresponding eigenvalues $`\{_{AA^{}}\}`$. Note that the allowed $`\{_{AA^{}}\}`$ are very much constrained by the constraints Eq. (8) on the $`n_r`$ themselves; however, we will not use the details of these explicitly.

First of all, observe that the A hexagons in turn form a triangular lattice, while the links $`rr^{}`$ between such hexagons can also be viewed as the links of this “A-lattice”, $`AA^{}rr^{}`$. In a given subsector with fixed $`\{_{AA^{}}\}`$, there remains one $`Z_N`$ degree of freedom for each such link. It is convenient to work in the number basis and label these remaining link degrees of freedom by

$`𝒩_{AA^{}}n_rn_r^{(0)}=(n_r^{}n_r^{}^{(0)})𝒩_{A^{}A},`$ (15)
where $`\{n_r^{(0)}\}`$ is one particular instance: $`_{AA^{}}=n_r^{(0)}+n_r^{}^{(0)}`$ (and our convention is that $`rA`$ and $`r^{}A^{}`$—see Fig. 3). Thus, $`𝒩_{AA^{}}𝒩_{AA^{}}`$ are oriented fields on the links of the A-lattice. The subsector is now completely specified by the conditions

$`{\displaystyle \underset{A^{}A}{}}𝒩_{AA^{}}=0,`$ (16)
which are the neutrality constraints Eq. (8) for the A hexagons.

The action of the Hamiltonian Eq. (9) in this subsector is readily described in terms of the new variables. Thus, the transverse field ($`h`$) terms are diagonal in the new number variables, while the B and C hexagon ring exchanges simultaneously raise (or lower) the three oriented number fields circulating around the corresponding A-lattice triangular plackets. Writing the raising operator for a given link number variable $`𝒩_{AA^{}}`$ as $`e^{i\mathrm{\Xi }_{AA^{}}}`$, the resulting Hamiltonian is

$`\widehat{H}[]=K{\displaystyle \underset{\mathrm{}}{}}(e^{i\mathrm{\Xi }_{AA^{}}}e^{i\mathrm{\Xi }_{A^{}A^{\prime \prime }}}e^{i\mathrm{\Xi }_{A^{\prime \prime }A}}+\mathrm{h}.\mathrm{c}.)`$ (17)
$`{\displaystyle \underset{AA^{}}{}}(\mathrm{\Gamma }_{AA^{}}e^{i(2\pi /N)𝒩_{AA^{}}}+\mathrm{h}.\mathrm{c}.),`$
where

$`\mathrm{\Gamma }_{AA^{}}=he^{i(2\pi /N)n_r^{(0)}}(1+e^{i(2\pi /N)_{AA^{}}}).`$ (18)
Together with the constraints Eq. (16), this is precisely the conventional $`Z_N`$ lattice gauge theory defined on the triangular A-lattice but with link-dependent $`\mathrm{\Gamma }_{AA^{}}`$ specific for the particular subsector $`\{_{AA^{}}\}`$. We can now use the conventional wisdom to characterize each such subsector and in turn the full hexagonal ring exchange Hamiltonian with $`K_A=0`$.

When $`h0`$, all the different subsectors are degenerate. The lowest energy state in each such subsector has the energy of $`2K`$ per triangle and is an equal weight superposition of all possible configurations of $`𝒩_{AA^{}}`$ that satisfy the constraints Eq. (16). Nonzero $`h`$ eliminates this degeneracy and selects one particular subsector, namely with all $`_{AA^{}}=0`$, as containing the true ground state of the full Hamiltonian with $`K_A=0`$. Indeed, treating $`\mathrm{\Gamma }_{AA^{}}`$ perturbatively, the lowest energy in a given subsector is

$`E_{GS}[]{\displaystyle \underset{\mathrm{}}{}}2K{\displaystyle \underset{AA^{}}{}}{\displaystyle \frac{|\mathrm{\Gamma }_{AA^{}}|^2}{2K[1\mathrm{cos}(2\pi /N)]}},`$ (19)
where for simplicity we assumed that the system has no boundaries. It is now clear that for small nonzero $`h`$ the ground state of the full ring exchange Hamiltonian with $`K_A=0`$ is in the subsector with all $`_{AA^{}}=0`$. The subsectors that are closest in energy have the smallest number of nonzero $`_{AA^{}}`$ and can be characterized as having alternating $`_{A_0A^{}}=+1`$ and $`_{A_0A^{}}=1`$ values on the six links to a given hexagon $`A_0`$ (this subsector is obtained from the ground state subsector by applying the hexagon ring exchange around the hexagon $`A_0`$). The energy gap to these subsectors is $`6h^2/K`$. We see that we have a peculiar situation where a nonzero transverse field $`h`$ is needed to stabilize this $`Z_N`$ deconfined ground state; this is because we are competing here against the $`Z_N\times Z_N`$ deconfined state that is obtained for large $`K_A,K`$.

We are all set to discuss confinement of charges in the model with $`K_A=0`$. The above analysis was carried out in the uncharged sector but is readily extended to the charged sectors. First, consider placing a pair of opposite charges on two B hexagons: $`_{rB_1}n_r=+1`$ and $`_{rB_2}n_r=1`$. Proceeding exactly as before, one is led to consider different subsectors (of this charged sector) specified by $`\{_{AA^{}}\}`$. In each such subsector, the Hamiltonian has precisely the form Eq. (17) with the number variables satisfying precisely the constraints Eq. (16). All information about the two charges is encoded in the allowed configurations $`\{_{AA^{}}\}`$, and one can clearly see that $`_{AA^{}}0`$ at least on a string of A-lattice bonds connecting $`B_1`$ and $`B_2`$. From the earlier arguments, the energy cost of introducing two such charges is then proportional to the length of this string, i.e., such charges are confined with the string tension $`h^2/K`$.

Consider now placing a pair of opposite charges on two A hexagons $`A_1`$ and $`A_2`$: $`_{rA_1}n_r=+1`$ and $`_{rA_2}n_r=1`$. The analysis of the subsectors $`\{_{AA^{}}\}`$ will be somewhat different in this case. For each such subsector in this charged sector there corresponds a subsector in the uncharged sector having exactly the same $`\{_{AA^{}}\}`$. It is convenient to “measure” each charged subsector relative to the corresponding uncharged subsector. This is achieved by defining link variables $`𝒩_{AA^{}}`$ via Eq. (15) using an uncharged instance $`\{n_r^{(0)}\}`$ of $`\{_{AA^{}}\}`$ (i.e., $`_{rR}n_r^{(0)}=0`$ for each $`R`$ and $`_{AA^{}}=n_r^{(0)}+n_r^{}^{(0)}`$ for each $`AA^{}`$). In each subsector, the Hamiltonian again has the form Eq. (17) when written in these link variables, which now satisfy new constraints $`_{A^{}A_1}𝒩_{A_1A^{}}=+1`$ and $`_{A^{}A_2}𝒩_{A_2A^{}}=1`$. This corresponding precisely to introducing two charges in the corresponding A-lattice gauge theory. Clearly, for large enough $`Kh`$, these charges will be deconfined.

We now have essentially complete description of the partially deconfined phase. Thus, one can readily identify the $`Z_N`$ vortex excitations of the A-lattice gauge theory with $`Z_N`$ vortices on the B and C hexagons. These vortices will have usual statistical interactions with the deconfined charges on the A hexagons. Also, as should become clear by reviewing the above discussion, we can essentially account for the different subsectors $`\{_{AA^{}}\}`$ by saying that there are additional particle excitations living on the A hexagons obtained from the ground state by the action of the corresponding A hexagon ring exchanges. These new particles have a “mass” of $`6h^2/K`$ and have no statistical interaction with the other particles.

We can now consider what happens when we allow nonzero $`K_A`$. As discussed above, the ring exchanges around the A hexagons introduce mixing between the different subsectors. However, as long as $`K_A`$ is much smaller than the corresponding gap $`h^2/K`$, the partially deconfined phase survives and is characterized by the same particle description.

Once $`K_A`$ is sufficiently large, the system enters the fully deconfined phase described in the main text.

## Appendix B Dual global $`Z_N`$ spin model

Here we summarize dual perspective on the hexagonal lattice ring exchange Hamiltonian Eq. (9). We work directly in the Hamiltonian language. Simple counting shows that the dimensionality of the physical Hilbert space is consistent with having one $`Z_N`$ degree of freedom per hexagon. Let us define

$`T_R^{}\psi _1^{}\psi _2\psi _3^{}\psi _4\psi _5^{}\psi _6,T_R^+(T_R^{})^{},`$ (20)
where we use the same sign convention as in Fig. 3. Let us also define

$`V_R^{}{\displaystyle \underset{R}{}}P^+,`$ (21)
where the product is along the vertical path that reaches $`R`$ as in Fig. 3. Note that the path “steps” through the same sublattice hexagons. \[If we were to take some other such path, we would need to replace some $`P^+`$ with $`P^{}`$. The total product is path-independent due to constraints Eq. (8).\] $`V_R^{}`$ can be thought of as a vortex creation operator.

We now interpret $`V_R^{}`$ as a $`Z_N`$ spin variable. It is easy to verify that $`T_R^+`$ is the corresponding conjugate variable (i.e., raising operatorZN ):

$`V_R^{}T_R^+=e^{i2\pi /N}T_R^+V_R^{}.`$ (22)
Also, we can readily “solve” for $`P_r^+`$:

$`P_r^+=V_{R_1}^{}V_{R_2}^{}V_{R_3}^{}.`$ (23)
The dual Hamiltonian is

$`H=K{\displaystyle \underset{R}{}}(T_R^++\mathrm{h}.\mathrm{c}.)h{\displaystyle \underset{\mathrm{}}{}}(V_{R_1}^{}V_{R_2}^{}V_{R_3}^{}+\mathrm{h}.\mathrm{c}.),`$ (24)
which is a global $`Z_N`$ spin model with three-spin interactions. A little thought shows that the model has in fact a $`Z_N\times Z_N`$ global symmetry corresponding to independent global rotations of the spins on two of the three sublattices. Note also that the three-spin interaction around triangles promotes ordering of the spins on the same sublattice. This is because two neighboring sites $`A`$ and $`A^{}`$ on the same sublattice share a $`BC`$ side in the respective triangle interactions $`\mathrm{}ABC`$ and $`\mathrm{}A^{}BC`$.

The global model clearly has a fully disordered phase for $`Kh`$. In the original ring exchange Hamiltonian, this corresponds to all vortices being gapped, and we obtain the $`Z_N\times Z_N`$ fully deconfined phase. Varying the A hexagon ring exchange coupling $`K_A`$ independently, for sufficiently small $`K_A`$ and large $`K`$ the system can clearly order on the A sublattice (i.e., A vortices condense), but remain disordered on the B and C sublattices. This is our partially deconfined $`Z_N`$ phase.


format: text

text:
# Heterogeneous decay of metastable phase on various centers -2

## 1 Intermediate approximate solution for several sorts of centers

In the case of nucleation on one sort of centers the expression for $`g_i(z)`$ is unknown only when the spectrum is cut off by exhaustion of heterogeneous centers. So, one doesn’t know it in the situation when the converging force of the operator $`S_i`$ is very strong. In the case of the several sorts of centers the situation is another. One doesn’t know every term in $`_jg_j(z)`$. Then one can come to the situation where according to the first iteration the spectrum is cut off by the exhaustion of supersaturation initiated by droplets formed on some sort of centers but in reality all heterogeneous centers of that sort are exhausted and those droplets practically don’t consume the vapor (due to their quantity). So, it is necessary to have some more precise expression for $`g_i`$ which allows to calculate the next iteration for $`\theta _i`$.

### 1.1 Mono-disperse approximation

The length corresponding to the cut-off of the spectrum by the substance exhaustion is practically one and the same for all sorts of the droplets (all $`\mathrm{\Gamma }_i`$ have approximately one and the same value111Alternative situations correspond to the asymptotic cases and have been already described.). Analyzing the subintergal expression in the equation for $`g_i`$ in the pseudo homogeneous case one can see that the subintegral expression $`s=(zx)^3\mathrm{exp}(\mathrm{\Gamma }_i_jg_j/\mathrm{\Phi }_{})`$ connected with a variation of the supersaturation in the formula for $`g_i`$ is a very sharp function of $`x`$. It is less than the function

$$s_{bel}=\mathrm{\Theta }(zx)(zx)^3$$
(1)
and greater than the function

$$s_{ab}=\mathrm{\Theta }(zx)(zx)^3\mathrm{exp}(\frac{\mathrm{\Gamma }_i_jf_j}{\mathrm{\Phi }_{}}\frac{x^4}{4}).$$
(2)
Fig. 6 illustrates the behavior of the values $`s_{ab}`$ and $`s_{bel}`$. It can be seen that they practically coincide. This is the reason of applicability of the first iteration as a good approximation for the precise solution for one sort of centers.

Let us introduce an approximation for this function. One has to extract the region of the sizes of the droplets which are essential in the vapor consumption. This consumption is essential when

$$z\mathrm{\Delta }x$$
(3)
where $`\mathrm{\Delta }x`$ is the length of the cut-off by the supersaturation. Certainly, the sizes of the droplets which are essential in the vapor consumption must be smaller than $`\mathrm{\Delta }x`$ because the effectiveness of the iteration procedure in the homogeneous case is based on the fact that the droplets formed at the almost ideal supersaturation determine the formation of the spectrum. One can see that $`s`$ is a very sharp function of $`x`$ and this is the reason for the monodisperce approximaion for $`s`$. Now one has to define characteristics of this approximation.

The condition for the differential width is that the subintegral function falls two times in comparison with an amplitude value. For the differential half-width $`\delta _{1/2}x`$ one can get the following expression

$$\delta _{1/2}x=(1\frac{1}{2^{1/3}})x.$$
(4)
The integral half-width $`\mathrm{\Delta }_{1/2}x`$ can be obtained from the corresponding equation

$$N_{ess}x^3=f_i\frac{x^4}{4}n_{\mathrm{}},$$
(5)
where $`N_{ess}`$ is the characteristic number of the droplets obtained as $`N_{ess}=f_i\mathrm{\Delta }_{1/2}xn_{\mathrm{}}`$ and the r.h.s. of the previous equation is the first iteration approximation. This gives

$$\mathrm{\Delta }_{1/2}x=\frac{1}{4}x$$
(6)
and it practically coincides with $`\delta _{1/2}x`$.

The subintegral function $`s`$ is now decomposed into the essential part where

$$x\frac{z}{4}$$

and the tail where

$$x\frac{z}{4}$$

The tail will be neglected and due to the small relative size of the essential part the mono-disperse approximation of sizes for the droplets in the essential part will be used. As the result one can get an approximation

$$g(z)=\frac{N(z/4)}{n_{\mathrm{}}}z^3,$$
(7)
where $`N(z/4)`$ is the number of the droplets appeared from $`x=0`$ till $`x=z/4`$.

Since the spectrum is cut off by the exhaustion of the supersaturation in a frontal (sharp) manner, the value of $`g_i`$ is unessential before $`z=\mathrm{\Delta }x`$ because it is small. After the moment of the cut-off it is unessential also because there is no formation of the droplets. So, instead of the previous approximation one can use

$$g_i(z)=\frac{N_i(\mathrm{\Delta }_ix/4)}{n_{\mathrm{}}}z^3.$$
(8)
Now let us turn to the heterogeneous case. The exhaustion of heterogeneous centers makes the subintegral function more sharp and the mono-disperse approximation becomes at $`x\mathrm{\Delta }_ix`$ even better than in the pseudo homogeneous situation. Certainly, one has to use $`N(\mathrm{\Delta }_ix/4)`$ calculated with account of exhaustion of heterogeneous centers (but at the coordinate, obtained without any account of exhaustion of heterogeneous centers).

Fig. 7 illustrates the form of the spectrum on the base of ”mono-disperse approximation”. The case of the pseudo homogeneous condensation is considered. Three curves are drawn: the spectrum in the mono-disperse approximation $`f_{appr}`$ ; the spectrum in the first iteration $`f_1`$ which can be considered as a very precise approximation and the numerical solution $`f`$. Two of them $`f`$ and $`f_1`$ coincide and correspond to a thick line. The back side in the first iteration is more sharp than in the ”mono-disperse approximation”. Nevertheless the deviation isn’t so essential. It can be eliminated by a perturbation theory.

One can get $`N(\mathrm{\Delta }_ix/4)`$ from the solution of the equations for the separate condensation process because only the lowest length of the cut-off is nesessary. This length is given without any cross influence taking into account due to the frontal character of the back side of the spectrum. Then one has to repeat this procedure with the already known $`N_i`$.

### 1.2 Final iterations

Now one can see a way to solve the system of the condensation equations.

A. Method of separate lengths.

At first one has to solve equations for all separate nucleation processes

$$g_i=f_i_0^z(zx)^3\mathrm{exp}(\mathrm{\Gamma }_i\frac{g_i}{\mathrm{\Phi }_{}})\theta _i𝑑xG_i(g_i,\theta _i),$$
(9)
$$\theta _i=\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}_0^z\mathrm{exp}(\mathrm{\Gamma }_i\frac{g_i}{\mathrm{\Phi }_{}})𝑑x)S_i(g_i)$$
(10)
for every $`i`$. This solution is given by the standard iteration procedure. In the second approximation one can calculate the value of $`\theta _{(2)}(\mathrm{\Delta }_ix/4)`$ explicitly

$$\theta _{i(2)}(\mathrm{\Delta }_ix/4)=\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}_0^{\mathrm{\Delta }_ix/4}\mathrm{exp}(\frac{\mathrm{\Gamma }_if_i}{4\mathrm{\Phi }_{}}z^4)𝑑x)=\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}(\frac{\mathrm{\Gamma }_if_i}{4\mathrm{\Phi }_{}})^{1/4}C),$$
(11)
where

$$C=_0^{1/4}\mathrm{exp}(z^4)𝑑x0.25.$$
(12)
Then

$$N_{i(2)}(\mathrm{\Delta }_ix/4)=\eta _{toti}(1\theta _{i(2)}(\mathrm{\Delta }_ix/4)).$$
(13)
One has to fulfil these calculations for every sort of the heterogeneous centers. Considering these approximations as the initial ones it is necessary to do only one more step of the iteration procedure to get the suitable results. They will be marked as the ’final’ ones. One has to calculate $`\theta _{ifinal}(\mathrm{})`$ and $`N_{finali}(\mathrm{})`$

$$\theta _{ifinal}(z)=\mathrm{exp}[f_i\frac{n_{\mathrm{}}}{\eta _{toti}}_0^z\mathrm{exp}(\frac{_j\mathrm{\Gamma }_iN_{j(2)}(\mathrm{\Delta }_jx/4)}{n_{\mathrm{}}\mathrm{\Phi }_{}}x^3)𝑑x],$$
(14)
$$N_{finali}(\mathrm{})=\eta _{toti}[1\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}(\frac{_j\mathrm{\Gamma }_iN_{j(2)}(\mathrm{\Delta }_jx/4)}{n_{\mathrm{}}\mathrm{\Phi }_{}})^{1/3}B].$$
(15)
These expressions are valid under the reasonable separation of heterogeneous centers into sorts when the centers with approximately one and the same height of the activation barrier are considered as one sort. Meanwhile it is obvious that when one splits one sort into very many sub-sorts then one can formally attain a wrong result.

Let us suppose that the number of heterogeneous centers ensures $`h_i1`$, i.e. the fall of the supersaturation leads to the interruption of formation of the droplets. One can split this sort into so many sub-sorts $`(j)`$ that for every sub-sort $`h_{ij}1`$, i.e. one can see the exhaustion of the sub-sort in the separate process of the condensation. Moreover, one can assume that this exhaustion is finished before $`z`$ attains the quarter of the length of the spectrum initiated by the fall of the supersaturation. Since $`\mathrm{\Delta }_{ij}x`$ is proportional to $`\eta _{ijtot}^{1/4}`$ (where $`\eta _{ijtot}`$ is the total number of the centers of the given sub-sort) and $`\delta _{ij}x`$ doesn’t depend on this quantity, the required property is obvious. Then after the summation over all sub-sorts one can see that the total number of the droplets formed up to $`\mathrm{\Delta }_ix/4`$ coincides with the total number of the heterogeneous centers. This conclusion is wrong.

The reason of the error is that the width of the spectrum is smaller than the width in the process of the separate formation. Evidently, the mono-disperse approximation doesn’t work at such distances. Under the reasonable definition of the heterogeneous sorts all characteristic lengths are different. So, ordinary, there is no such effect in this situation. But still one has to correct these results. The simplest way is to consider all centers with the difference in the activation barrier height less than one thermal unit as one and the same sort. This solves the problem, but this requirement is rather artificial.

Now these results will be corrected in more elegant way. The mono-disperse approximation leads to some already defined functional form for $`g`$ and for $`\zeta `$. The form is already known and now only parameters in these functional dependencies have to be determined.

B. Method of the unique length.

The unique length of the spectrum $`\mathrm{\Delta }x`$ will be chosen. At this very moment $`\mathrm{\Delta }x`$ is unknown, but satisfies the following inequality:

$$\mathrm{\Delta }x\mathrm{\Delta }_ix$$

for every sort of the heterogeneous centers. This leads to

$$g_i(z)=\frac{N_i(\mathrm{\Delta }x/4)}{n_{\mathrm{}}}z^3,$$
(16)
which gives

$$N_{i(2)}(\mathrm{\Delta }x/4)=\eta _{itot}(1\theta _{i(2)}(\mathrm{\Delta }x/4)),$$
(17)
$$\theta _{i(2)}(\mathrm{\Delta }x/4)=\mathrm{exp}[f_i\frac{n_{\mathrm{}}}{\eta _{toti}}_0^z\mathrm{exp}(\frac{\mathrm{\Gamma }x^3}{\mathrm{\Phi }_{}n_{\mathrm{}}}\underset{j}{}N_j(\mathrm{\Delta }x/4))𝑑x].$$
(18)
The weak dependence $`\mathrm{\Gamma }_i\mathrm{\Gamma }`$ on the sort of the centers is neglected here for simplicity, while the strong dependence $`f_i`$ on the sort of centers is taken into account. After the substitution one can get the system of the algebraic equations for $`N_i(\mathrm{\Delta }x/4)`$

$$N_i(\mathrm{\Delta }x/4)=\eta _{itot}(1\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}_0^{\mathrm{\Delta }x/4}\mathrm{exp}(\frac{\mathrm{\Gamma }}{\mathrm{\Phi }_{}}\frac{x^3}{n_{\mathrm{}}}\underset{j}{}N_j(\mathrm{\Delta }x/4))𝑑x)).$$
(19)
Let us simplify the last system. To calculate the integral one can note that

$$_0^x\mathrm{exp}(x^3)𝑑xxx1/4.$$
(20)
So, the trivial dependence of the r.h.s. on $`N_j`$ disappears. One can come to

$$N_i(\mathrm{\Delta }x/4)=\eta _{itot}(1\mathrm{exp}(f_i\frac{n_{\mathrm{}}}{\eta _{toti}}\frac{\mathrm{\Delta }x}{4})).$$
(21)
On the other hand one can use the sense of $`\mathrm{\Delta }x`$ as the half-width of the spectrum due to the fall of the supersaturation:

$$\frac{\mathrm{\Delta }x^3}{n_{\mathrm{}}}\underset{j}{}N_i(\mathrm{\Delta }x/4)\frac{\mathrm{\Gamma }}{\mathrm{\Phi }_{}}=1.$$
(22)
After the substitution one can get the equation for $`\mathrm{\Delta }x`$

$$\frac{\mathrm{\Delta }x^3}{n_{\mathrm{}}}\underset{j}{}\eta _{jtot}(1\mathrm{exp}(f_j\frac{n_{\mathrm{}}}{\eta _{totj}}\frac{\mathrm{\Delta }x}{4}))\frac{\mathrm{\Gamma }}{\mathrm{\Phi }_{}}=1.$$
(23)
To solve this algebraic equation one can choose the initial approximation as a solution of the following equation

$$\frac{\mathrm{\Delta }x^4}{4}\underset{j}{}f_j\frac{\mathrm{\Gamma }}{\mathrm{\Phi }_{}}=1.$$
(24)
Iterations are constructed according to

$$\frac{\mathrm{\Delta }_{(i)}x^3}{n_{\mathrm{}}}\underset{j}{}\eta _{jtot}(1\mathrm{exp}(f_j\frac{n_{\mathrm{}}}{\eta _{totj}}\frac{\mathrm{\Delta }_{(i1)}x}{4}))\frac{\mathrm{\Gamma }}{\mathrm{\Phi }_{}}=1.$$
(25)
When the last equation is solved then equation (21) expresses $`N_i(\mathrm{\Delta }x/4)`$ through $`\mathrm{\Delta }x`$ and solves the problem. Equation (16) gives the expression for $`g_i`$ and for the supersaturation as the function of time.

One can estimate the relative error by the pseudo homogeneous case which leads to

$$\frac{N_i(\mathrm{})N_{finali}(\mathrm{})}{N_i(\mathrm{})}\frac{AB}{A}0.02$$
(26)
The investigation of the process of nucleation on several sorts of heterogeneous centers is completed.

## 2 Solution for heterogeneous centers with quasi continuous activity

### 2.1 Iteration procedure

The smoothness of $`\eta (w)`$ is supposed here. It allows to hope that the cross influence of exhaustion of one type of heterogeneous centers on the cut-off of the droplets formation on another sort of centers will be not so strong and one can apply some modifications of the standard iteration procedure. It can be constructed in a following way: for initial approximations one can choose

$$g_0(z,w)=0,\theta _0=1;$$
(27)
the recurrent procedure is defined according to

$$g_{i+1}(z,w)=G_w(g_i^{tot},\theta _i(w)),g_i^{tot}(z)=dwg_i(z,w).$$
(28)
$$\theta _{i+1}(z,w)=S_w(g_i^{tot})$$
(29)
The monotonic properties of $`G_w`$, $`𝑑w`$ and $`S_w`$ lead to the chains of inequalities analogous to (LABEL:chain)-(LABEL:chain), which ensure estimates for precision of approximations. Certainly, to calculate iterations one has to use some expression for $`\eta _{tot}(w)`$. So, at first one has to consider the limit cases and obtain some estimates.

One can see the absence of $`\eta _{tot}`$ in the r.h.s. of (29). So, for arbitrary $`\eta _{tot}`$ the power of exhaustion will be determined only by $`w`$. The centers with the high activity $`w`$ are almost exhausted during the process of the condensation. The centers with the relatively small activity $`w`$ remain free. The intermediate region has rather small size. To give the qualitative estimates one can see that the supersaturation $`\zeta `$ appears only in the function $`f_\zeta `$. For $`f_\zeta `$ one can obtain the following estimate

$$f_\zeta f_{}\mathrm{exp}(constx^ϵ)3ϵ4$$
(30)
This estimate goes from the obvious fact that the spectrum of the droplets must be wider than the mono-disperse one and the intensity of nucleation must decrease in time. On the base of the last estimates one can see that the width $`\mathrm{\Delta }x`$ of the spectrum is a well defined value.

The rough approximation for exhaustion of heterogeneous centers can be obtained from the first approximation in the iteration procedure

$$\theta =\mathrm{exp}[f_{}\mathrm{exp}(\lambda w)n_{\mathrm{}}z]$$
(31)
and for the value at the end of nucleation period

$$\theta _{final}=\mathrm{exp}[f_{}\mathrm{exp}(\lambda w)n_{\mathrm{}}\mathrm{\Delta }x].$$
(32)
One can define $`w_<`$ and $`w_>`$ according to

$$w_<=w_0\frac{ϵ}{\lambda },w_>=w_0+\frac{ϵ}{\lambda },$$
(33)
where

$$w_0=\frac{1}{\lambda }\mathrm{ln}(\frac{\mathrm{ln}2}{f_{}n_{\mathrm{}}\mathrm{\Delta }x}),ϵ1.$$
(34)
These constructions are illustrated in Fig. 8.

For $`w>w_>`$ almost all heterogeneous centers are exhausted. For $`w<w_<`$ almost all heterogeneous centers remain free. This remark allows to rewrite the expression for $`g^{tot}`$ in the following form

$$g^{tot}=_{w_<}^{w_>}𝑑wg(z,w)+_{w_>}^{\mathrm{}}\eta _{tot}(w)𝑑w\frac{z^3}{n_{\mathrm{}}}.$$
(35)
The size of the intermediate region corresponds to a variation of $`w`$ of an order $`1/\lambda `$ or to the variation of $`\mathrm{\Delta }F`$ in one thermal unit. Since $`\mathrm{\Delta }F1`$, the relative variation of $`\mathrm{\Delta }F`$ in the intermediate region is small. So, it is reasonable to put in this region

$$\eta _{tot}(w)=\eta _{}=const.$$
(36)
One can spread this approximation over the region $`w>w_>`$ because the behavior of $`\eta _{tot}`$ for these centers isn’t important. Only $`\eta _{tot}𝑑w`$ over this region is essential. Certainly, one has to get the boundary $`w_{max}`$ of this region by

$$_{w_>}^{w_{max}}\eta _{}𝑑w=\eta _{}(w_{max}w_>)=_{w_>}^{\mathrm{}}\eta _{tot}(w)𝑑w.$$
(37)
For $`w<w_<`$ the accuracy of this approximation is not so essential because these centers remain free.

After the scale transformations $`\mathrm{\Gamma }g^{tot}/\mathrm{\Phi }_{}G;\mathrm{\Gamma }g/\mathrm{\Phi }_{}g`$

$$g=\frac{\mathrm{\Gamma }f_{}\eta _{}}{\mathrm{\Phi }_{}}\mathrm{exp}(\lambda w)_0^z(zx)^3\mathrm{exp}(G)\theta 𝑑x,$$
(38)
$$G=_{\mathrm{}}^{w_{max}}𝑑wg,$$
(39)
$$\theta =\mathrm{exp}[f_{}n_{\mathrm{}}\mathrm{exp}(\lambda w)_0^z\mathrm{exp}(G)𝑑x].$$
(40)
By the appropriate shift of $`w`$ one can put $`w_{max}`$ to the zero. In the appropriate scale of $`w`$ one can put $`\lambda =1`$. In the appropriate scale of $`z`$ one puts the coefficient in (38) to unity. As the result, only one parameter - the coefficient $`f_{}^{3/4}\mathrm{exp}(3\lambda w_{max}/4)\mathrm{\Phi }_{}^{1/4}/(\mathrm{\Gamma }^{1/4}\eta _{}^{1/4})`$ remains in the last equation. One can mark it by $`A`$ and get

$$g=\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G)\theta 𝑑x,$$
(41)
$$G=_{\mathrm{}}^0𝑑wg,$$
(42)
$$\theta =\mathrm{exp}[A\mathrm{exp}(w)_0^z\mathrm{exp}(G)𝑑x].$$
(43)
An iteration procedure can be constructed as

$$g_{i+1}=\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G_i)\theta _i𝑑x,$$
(44)
$$G_i=_{\mathrm{}}^0𝑑wg_i,$$
(45)
$$\theta _{i+1}=\mathrm{exp}(A\mathrm{exp}(w)_0^z\mathrm{exp}(G_i)𝑑x),$$
(46)
$$g_0=0,\theta _0=1.$$
(47)
In the first iteration

$$g_1=\mathrm{exp}(w)\frac{z^4}{4},$$
(48)
$$G_1=_{\mathrm{}}^0𝑑wg_1=\frac{z^4}{4},$$
(49)
$$\theta _1=\mathrm{exp}(A\mathrm{exp}(w)z).$$
(50)
The second approximation gives for $`\theta `$ the following result

$$\theta _2=\mathrm{exp}[A\mathrm{exp}(w)_0^z\mathrm{exp}(\frac{z^4}{4})𝑑z],$$
(51)
then for the final value

$$\theta _{final2}=\mathrm{exp}[A\mathrm{exp}(w)1.28^{1/4}]$$
(52)
and for $`N^{tot}`$

$$N^{tot}=_{\mathrm{}}^0[1\mathrm{exp}(Aexp(w)1.28^{1/4})]𝑑w\eta _{}.$$
(53)
For $`g`$ in the second approximation

$$g_2=\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(A\mathrm{exp}(w)x)\mathrm{exp}(x^4/4)𝑑x.$$
(54)
The analytical expression for $`g_2`$ and further iterations can not be calculated.

### 2.2 Universal solution

The system of equations (41) - (43) doesn’t allow the universal solution as in the case of homogeneous condensation when (41) hasn’t $`\theta `$ in the r.s.h.. But the operator in the r.h.s. of (41) ensures a rapid convergence of the iterations. The worst situation for the iteration convergence is $`A=0`$. In this situation one has after rescaling the universal system

$$g=\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G)𝑑x,$$
(55)
$$G=_{\mathrm{}}^0g𝑑w.$$
(56)
This system has the universal solution222The trivial integration over $`w`$ leads to the equation analogous to the homogeneous case. which will be marked by $`G_0`$. This universal form of the size spectrum is shown in Fig. 9.

The evolution is determined by the first three momentums (and the zero one) of the distribution function

$$\mu _i(z)=_0^zx^i\mathrm{exp}(G)\theta dx,i=0,1,2,3$$
(57)
and in the pseudo homogeneous case it is determined by

$$\mu _i(z)=_0^zx^i\mathrm{exp}(G_0)𝑑x.$$
(58)
After the end of the short period of intensive formation of droplets one can substitute in (57) and (58) $`\mathrm{}`$ instead of $`z`$ in the region of the integration. So, the further evolution will be determined by four constants333 The main influence is due to $`\mu _0`$. $`\mu _i(\mathrm{})i=0,1,2,3`$. If one takes for $`G`$ the universal solution $`G_0`$ (when $`A=0`$) then the values of $`\mu _i`$ are the universal constants which can be obtained by the unique numerical solution of the last system.

Now one has to return to the iteration procedure. One can use $`G_0`$ as the initial approximation in the iteration procedure and get

$$\theta _1=\mathrm{exp}[A\mathrm{exp}(w)_0^z\mathrm{exp}(G_0)𝑑x],$$
(59)
$$g_2=\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G_0)\mathrm{exp}(A\mathrm{exp}(w)_0^x\mathrm{exp}(G_0)𝑑x^{})𝑑x,$$
(60)
$$G_2=_{\mathrm{}}^0\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G_0)\mathrm{exp}(A\mathrm{exp}(w)_0^x\mathrm{exp}(G_0)𝑑x^{})𝑑x𝑑w.$$
(61)
The decomposition of exponents gives

$$G_2=\underset{i=0}{\overset{\mathrm{}}{}}\frac{(A)^i}{(i+1)!}P_i(z),$$
(62)
where

$$P_i=_0^z(zx)^3\mathrm{exp}(G_0)J_{00}^i(x)𝑑x,$$
(63)
$$J_{00}(x)=_0^x\mathrm{exp}(G_0)𝑑x.$$
(64)
Then the function $`B(z)=_0^z\mathrm{exp}(G_2(x))𝑑x`$ at $`z=\mathrm{}`$ will be the following

$$B(z=\mathrm{})=A_0^{\mathrm{}}\mathrm{exp}(G_2)𝑑x=A_0^{\mathrm{}}\mathrm{exp}(P_0(z))\mathrm{exp}(\underset{i=1}{\overset{\mathrm{}}{}}\frac{(A)^i}{(i+1)!}P_i(z))𝑑z.$$
(65)
The decomposition of the last exponent leads to

$$B(z=\mathrm{})=\underset{i=1}{\overset{\mathrm{}}{}}\underset{j=0}{\overset{\mathrm{}}{}}\frac{(1)^j}{j!}A\frac{(A)^{ij}}{((i+1)!)^j}C_{ij},$$
(66)
where

$$C_{ij}=_0^{\mathrm{}}\mathrm{exp}(P_0(x))P_i^j𝑑x$$
(67)
are the universal constants. For the total number of the droplets with a given activity

$$N_{tot}(w)=\eta _{}(1\mathrm{exp}(B(z=\mathrm{})\mathrm{exp}(w))).$$
(68)
In the same manner one can take into account the deviation of $`\eta _{}`$ from a constant value. After the decomposition of $`\eta _{}`$ into Tailor series one can take the initial approximation as $`G_0`$ and act in the manner similar to (59) - (68). As a result, some expansions in the powers of the parameter $`A`$, the derivatives of $`\eta _{tot}`$ and the universal constants can be seen.

Fig. 10 illustrates these constructions.

Here the most interesting function $`I_{00}=_0^{\mathrm{}}\mathrm{exp}(G)𝑑x`$ is drawn. Two curves are shown. the lower curve corresponds to the precise numerical solution, i.e. to $`_0^{\mathrm{}}\mathrm{exp}(G)𝑑x`$. The upper curve corresponds to $`B(z=\mathrm{})/A=_0^z\mathrm{exp}(G_2)𝑑x`$. This curve is very close to numerical solution. It is very important that both curves resemble a straight line. So, already the linear approximation on $`A`$ corresponding to the account of the first leading (linear) term will be good. It means that in the last expression one can take instead of $`B(z=\mathrm{})`$ the value $`A_0^{\mathrm{}}\mathrm{exp}(G_0)𝑑x`$ (on the base of the first iteration one puts $`A_0^{\mathrm{}}f_1𝑑x`$).

### 2.3 Wide spectrum of active centers

Three regions were extracted in the spectrum of activities: the region of active centers (they are almost exhausted), the region of centers with small activity (they remain practically free during the whole process) and the intermediate region. In the majority of the situations the intermediate region has relatively small size in comparison with the active region (the passive region has no size - it can be spread to infinity). Then one can fulfill some further simplifications.

One can transmit the point $`w=0`$ to the special activity for which

$$\theta _{w=0}(\mathrm{})=\frac{1}{2},$$

i.e. the half of the centers became the centers of droplets. Then the system of the condensation equations after the rescaling of $`w`$ can be written as

$$g=a_0\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G)\theta 𝑑x,$$
(69)
$$G=_{\mathrm{}}^{w_+}g(w)𝑑w,$$
(70)
$$\theta =\mathrm{exp}(a_1\mathrm{exp}(w)_0^z\mathrm{exp}(G)𝑑x),$$
(71)
where $`w_+`$ is the upper boundary of spectrum and $`a_0,a_1`$ are some constants. Having rescaled $`z,x`$ one can put $`a_0=1`$. The condition of the choice for $`w=0`$ gives:

$$a_1=\frac{\mathrm{ln}2}{_0^{\mathrm{}}\mathrm{exp}(G)𝑑x}.$$
(72)
So, there remains only one parameter $`w_+`$.

Let us explicitly extract the intermediate and the active regions. One can see that for two different activities $`w_1`$ and $`w_2`$ the following equation is valid

$$\frac{\mathrm{ln}(\theta (w_1,x))}{\mathrm{exp}(w_1)}=\frac{\mathrm{ln}(\theta (w_2,x))}{\mathrm{exp}(w_2)}.$$
(73)
The same is valid for the final values of $`\theta `$. The value of $`\mathrm{ln}(\theta (w,x))/\mathrm{exp}(w)`$ is invariant for the different activities in one and the same process. So, one can put $`w=w_{++}2÷3`$ as the boundary between the active region and the intermediate region. In the same manner one can separate the region of the passive centers by the boundary $`w_{}=w_{++}`$. One can neglect the substance in the droplets on the passive centers and get444The absence of the active and intermediate regions ($`w_+<w_{}`$) means that the condensation occurs in the pseudo homogeneous way.

$$G_w_{}^{w_+}g(w)𝑑w.$$
(74)
One can denote by $`n_{\mathrm{}}G_+`$ the number of the substance molecules in the droplets formed on the active centers. The system of the condensation equations can be rewritten as

$$g=a_0\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G)𝑑x,$$

$$G=_{\mathrm{}}^{w_{++}}g(w)𝑑w+G_+,$$

$$\theta =\mathrm{exp}(a_1\mathrm{exp}(w)_0^z\mathrm{exp}(G)𝑑x).$$

Now one can write an evident expression for $`G_+`$

$$G_+=\frac{\eta _{}}{n_{\mathrm{}}}(w_+w_{++})z^3.$$
(75)
The system of the condensation equations now is the following one

$$g=a_0\mathrm{exp}(w)_0^z(zx)^3\mathrm{exp}(G)𝑑x,$$

$$G=_{\mathrm{}}^{w_{++}}g(w)𝑑w+\frac{\eta _{}}{n_{\mathrm{}}}(w_+w_{++})z^3,$$

$$\theta =\mathrm{exp}(a_1\mathrm{exp}(w)_0^z\mathrm{exp}(G)𝑑x).$$

The value of $`w_{++}`$ is universal, but the coefficient in the term $`\eta _{}(w_+w_{++})z^3/n_{\mathrm{}}`$ depends on parameters.

Having defined

$$G_{}=_{\mathrm{}}^{w_{++}}g(w)𝑑w,$$
(76)
one can propose the following estimate

$$G_{}(w_{++}w_{})\frac{\eta _{}}{n_{\mathrm{}}}z^3.$$
(77)
When

$$w_+w_{++}w_{++}w_{}4÷6$$

one can approximately get

$$G=\frac{\eta _{}}{n_{\mathrm{}}}(w_+w_{++})z^3.$$

Then the first equation isn’t necessary at all. So, there is no necessity to put $`a_0`$ to $`1`$ and one can use the arbitrary scale to cancel another constant. The system looks like

$$\theta =\mathrm{exp}(a_1\mathrm{exp}(w)_0^z\mathrm{exp}(G)𝑑x),$$

$$G=\frac{\eta _{}}{n_{\mathrm{}}}(w_+w_{++})z^3.$$

Having rescaled $`x,z`$, one can put $`\eta _{}(w_+w_{++})=n_{\mathrm{}}`$ and the system looks like the universal expression:

$$\theta =\mathrm{exp}(\frac{\mathrm{ln}2}{_0^{\mathrm{}}\mathrm{exp}(z^3)𝑑x}\mathrm{exp}(w)_0^z\mathrm{exp}(z^3)𝑑x).$$

The behavior of $`\theta `$ as a function of $`z`$ for the different values of $`w`$ is shown in Fig. 11.

One has to stress that in the general case of quasi continuous spectrum all results concerning the intermediate approximate solution can be used. One has only to substitute the sum over all sorts of centers by the integral over activities with appropriate normalizing coefficient: $`_idw`$.

Moreover, the procedure of intermediate approximate solution can be directly applied to the situation with a mixed quasi continuous spectrum of activities and some discrete $`\delta `$-like peaks for several sorts of heterogeneous centers.

This remark completes the description of the nucleation for the decay on heterogeneous centers with different activities. Unfortunately there are no experimental investigations which can be compared with the theory. It is rather natural because earlier there were no formulas which could give any physical information and experimental investigations would give only some chaotic experimental data. Now one can easily check this theory. This publication can be regarded as an invitation to open experimental investigations in this field. One can use the same devices as for the ordinary homogeneous nucleation investigations and calculate the number of droplets. Experiments on this theory can give information about the properties of heterogeneous centers in complex systems. The special interest is to solve some extreme problems to produce mainly the droplets on a given sort of centers when there are several sorts of centers in the system. With the help of this theory this prolem can be solved.


text:
# Expanded version of the contribution in NEUTRINO2002, May 25-30,2002, Munich, Germany Neutralino relic density from direct non-equilibrium production and intermediate scales

## 1 Introduction

In a series of works , various implications for an intermediate unification scale in supersymmetric models have been considered. In particular, it was pointed out that the lightest neutralinos (LSP)–nucleon cross sections are very sensitive to the variation of this scale. For instance, for $`M_I10^{12}10^{14}`$ GeV the cross sections are enhanced and a large region of the parameter space of the minimal supersymmetric standard model (MSSM) is compatible with the sensitivity of the current and near future direct detection dark matter experiments (DAMA, CDMS, EDELWEISS). In fact, such results can also be obtained within the GUT unification scenario, i.e., $`M_I10^{16}`$ GeV, but it requires tuning the available parameter space: either large $`\mathrm{tan}\beta `$ ($`\mathrm{tan}\beta 25`$) or a specific choice for non–universal soft SUSY breaking terms so that the value of $`\mu `$ term is reduced and hence a larger Higgsino components of the LSP is obtained, which is essential for enhancing the LSP–nucleon cross sections .

One problem is than in GUT unification scenarios with standard Big Bang (BB) cosmological expansion where the LSP is produced in a thermal environment dominated by radiation, such large cross section leads to a very small LSP relic density (outside the cosmological and astrophysical interesting region $`0.1\mathrm{\Omega }h^20.3`$). In this scheme there is a generic anti-correlation between the LSP relic density and its quark or nucleon cross sections: one usually obtains

$$\mathrm{\Omega }_\chi h^2\frac{C}{\sigma _\chi ^{\mathrm{ann}}v},$$
(1)
where $`\sigma ^{\mathrm{ann}}`$ is the annihilation cross section for neutralino pair, $`v`$ is the pair relative velocity, and $`..`$ denotes thermal averaging. Scattering ($`\sigma _{\chi q}`$) and annihilation $`\sigma _\chi ^{\mathrm{ann}}`$ LSP cross sections can be related by a cross symmetry. Therefore, a large scattering cross section $`\sigma _{\chi q}`$ leads generically to a large annihilation cross section and a small relic density as Eq. 1 suggest. Indeed, in these models, it was observed that LSP–proton cross sections of order $`10^6`$ Pb correspond to too low relic densities ($`<0.0050.01`$ ).

We want to address here the relic density question: in the context of diverse non-standard cosmological scenarios, what is the effect of changing the unification scale from GUT scale ($`10^{16}`$ GeV) down to intermediate scale ($`10^{12}`$ GeV) and whether we still have any anti-correlation between $`\mathrm{\Omega }_\chi h^2`$ and $`\sigma _\chi ^{\mathrm{ann}}`$.

## 2 The standard BB scenario

The LSP relic density is given by

$$\mathrm{\Omega }_\chi h^2=\rho _\chi /\rho _{crit},$$

where $`\rho _\chi `$ is the LSP mass density and $`\rho _{crit}=3H_0^2/8\pi G1.9\times 10^{29}h^2g/cm^3`$. The evolution of the number density $`n_\chi =\rho _\chi /m_\chi `$ can be described to a precision sufficient for our purposes by the Boltzmann equation

$$\frac{dn_\chi }{dt}+3Hn_\chi =\sigma _\chi ^{ann}v\left[(n_\chi )^2(n_\chi ^{eq})^2\right],$$
(2)
where $`H`$ is the Hubble expansion rate. In the standard calculation, it is assumed that in the early universe $`(Tm_\chi `$) the LSP density was initially equal to its equilibrium value $`n_\chi =n_\chi ^{eq}`$. As the universe expands the LSP density traces its equilibrium value. When the particle becomes non relativistic ( $`T<m_\chi `$) the number density of $`\chi `$’s drops exponentionally. The annihilation rate $`\mathrm{\Gamma }_\chi \sigma _\chi ^{ann}vn_\chi `$ becomes smaller than the expansion rate, $`\mathrm{\Gamma }_\chi <H`$. At this point, the LSP decouples and fall out from being in equilibrium with the radiation, their number remains constant after that and its density is only diluted by the universe expansion.

As long as the LSP is non-relativistic, we can expand the averaged cross section $`\sigma _\chi ^{ann}v`$ as follows

$$\sigma _\chi ^{ann}v=\alpha _s+\alpha _pv^2,$$
(3)
where $`\alpha _s`$ is the s-wave contribution at zero relative velocity and $`\alpha _p`$ contains contribution from both s- and p- waves. The relic density solution to the Boltzamn equation is accurately given by the approximate formula (c.f. Eq.1)

$$\mathrm{\Omega }_\chi h^2\frac{8.8\times 10^{11}GeV^2}{\sqrt{g_{}}(\alpha _s/x_F+3\alpha _p/x_F^2)},$$
(4)
where the freezing-out temperature $`x_Fm_\chi /T_F`$ can be estimated as

$$x_F\mathrm{ln}\frac{c(\alpha _s+6\alpha _p/x_F)}{\sqrt{g_{}x_F}}.$$
(5)
In these formulas $`g_{}`$ is the effective number of relativistic degrees of freedom at $`T_F`$ and $`c`$ is constant of order $`10^1`$ .

## 3 Alternative scenarios

In alternative cosmological scenarios, it is often assumed the existence of one (o more) epoch before the radiation-dominated era of the universe where the energy is dominated by coherent oscillating fields, the so called reheating era. A reheating process is associated by example with the final stage of Inflation. It can be also, independently or not, associated with the oscillations of moduli or Polony fields appearing in string theory after flat directions acquire a mass from SUSY breaking.

The details and physical origin of these reheating are largely unimportant. Its effects can be characterized by a reheating temperature $`T_{RH}`$ which is conventionally defined (see the definition below) as the temperature at which the oscillating field energy ceases to dominate the cosmological evolution and starts the radiation dominated epoch.

The standard BB scenario presented in the previous section can be seen as an special case when $`T_{RH}`$ is too much higher than $`T_F`$ (of order $`10^{910}`$ GeV ). In this case, the reheating epoch has no relevance in the final output of the relic density.

A problem could appear if the decay of the oscillating field occurs after or during the primordial BBN nucleosynthesis. This is drastically solved imposing from the onset that the reheating temperature is larger than $`1`$ MeV.

Finally, in interesting intermediate cases, a relatively low value of $`T_{RH}`$ such that $`T_F>T_{RH}`$ can have qualitative and quantitative implications on the predictions of the relic abundance of the LSP as discussed in Ref. . The main reason for this is that now the decoupling of the LSP particle occurs in an expanding universe dominated by matter instead than radiation and that has effects. We will see below that theories with intermediate unification scale predict low reheating temperature ($`<GeV`$ and below) and they can be considered as interesting scenarios.

It is useful to use the following expression which can be interpreted as a quantitative definition of the temperature $`T_{RH}`$ and which connects it to the decay width of the oscillating field:

$$\mathrm{\Gamma }_\varphi =\left(\frac{4\pi ^3g_{}(T_{RH})}{45}\right)^{1/2}\frac{T_{RH}^2}{M_P}.$$
(6)
Given any concrete theory, the reheating temperature can be related to the their physical parameters inverting Eq.6 for $`T_{RH}`$ as we can see on the following.

Let us consider first an inflation scenario. The inflaton $`(\varphi )`$ decay width which $`\mathrm{\Gamma }_\varphi `$ can be written, in some some supersymmetric hybrid inflation which has been considered , as

$$\mathrm{\Gamma }_\varphi =\frac{1}{8\pi }\left(\frac{M_f}{\varphi }\right)^2m_\varphi ,$$

where $`m_\varphi `$ is the inflaton mass. Here $`M_f`$ is the mass of the particle $`f`$ that the inflaton decay to (i.e. a right handed neutrino or a sneutrino). $`M_f`$ should be less than the inflaton mass to allow for the decay $`\varphi ff`$. Finally $`\varphi `$ is the vev of the inflaton field which is of the order of the unification scale.

In Ref.  it has been shown that a intermediate unification scale of order $`M_I10^{1112}`$ GeV is favored by inflation. In this case, recalling that the inflaton mass is constrained by

$$m_\varphi <M_I^2/M_P,$$

one obtains ($`\varphi 1.7\times 10^{16}`$ GeV ) an inflaton mass of the order $`m_\varphi 10^2`$ GeV. For the ’standard GUT scenario’ with $`M_I10^{16}`$ one obtains however $`m_\varphi 10^{13}`$ GeV . From Eq.(6) we find that $`T_{RH}𝒪(1)`$ GeV; a value much smaller than the values quoted in the standard GUT scenario ($`T_{RH}10^{11}`$). This reheating temperature is of the same order or lower than the typical freeze-out temperature $`T_Fm_\chi /20O(15)`$ GeV.

A detailed analysis of the relic density with a low reheating temperature has been considered in Ref.  with different possible generic scenarios. The first scenario corresponds to a case where the LSP is never in chemical equilibrium either before or after reheating. In the second the LSP’s reach chemical equilibrium but they freeze out before the completion of the reheat process $`T_F>T_{RH}`$. Each of these scenarios leads to different qualitative and quantitative predictions for the relic density which can become quite different from the standard computation we summarized in Eq.(4).

CASE A) In case of non–equilibrium production and freeze out at the early times and supposing that the number density of the LSP $`n_\chi `$ is much smaller that $`n_\chi ^{eq}`$, the Boltzmann equation (2) can be approximated and solved. In this case one gets

$`\mathrm{\Omega }h^2=2.1\times 10^4\lambda _{}\left(\alpha _s+{\displaystyle \frac{\alpha _p}{4}}\right),`$ (7)
where

$`\lambda _{}={\displaystyle \frac{g^2}{4}}\left({\displaystyle \frac{g_{}(T_{RH})}{10}}\right)^{3/2}\left({\displaystyle \frac{10}{g_{}(T_{})}}\right)^3{\displaystyle \frac{(10^3T_{RH})^7}{m_\chi ^5}},`$ (8)
and where $`g`$ is the number of degrees of freedom of the LSP and $`T_{}`$ is the temperature at which most of the LSP production takes place, it is given by $`T_{}4m_\chi /15`$. As we can see this $`\mathrm{\Omega }h^2`$ is proportional to the annihilation cross section, instead of being inversely proportional, as in Eq.(4). However, the assumption that $`n_\chi <<n_\chi ^{eq}`$ leads to a sever constraint on the annihilation cross section : $`\alpha _s<\overline{\alpha _s}`$ and $`\alpha _p<\overline{\alpha _p}`$ where $`\alpha _s10^{10}10^9`$ and $`\alpha _p10^910^8`$. So that for large neutralino cross sections, which we are interested in, eq.(7) can not be applied.

CASE B) With large annihilation cross section the LSP reaches equilibrium before reheating as discussed in Ref. and its relic density $`\mathrm{\Omega }h^2`$ is very close to the one obtained by using the standard calculation of relic density in Eq.(4). In this case, $`T_F`$ is obtained by solving a different equation as before

$`x_F=\mathrm{ln}\left[C\left(\alpha _sx_F^{5/2}+{\displaystyle \frac{5}{4}}\alpha _px_F^{3/2}\right)\right],`$ (9)
where $`C`$ is an unimportant constant now. The relic density $`\mathrm{\Omega }h^2`$ is given by

$$\mathrm{\Omega }h^2=\frac{g_{}^{1/2}(T_{RH})}{g_{}(T_F)}\frac{2.3\times 10^{11}T_{RH}^3GeV^2}{m_\chi ^3(\alpha _sx_F^4+4\alpha _px_F^5/5)}.$$
(10)
In this case the relic density is again inversely proportional to the annihilation cross section as in Eq.(4). Moreover, it has a further suppression due to the very low reheating temperature effect. The typical value for the relic density in this case for $`\sigma ^{ann}10^6`$ GeV-2 is $`\mathrm{\Omega }h^2<10^3`$. So it is even worse result than the BB standard computation’s one.

## 4 Our Model

In superstring theory, which is the first motivation for the idea of considering an intermediate scale , there are many moduli fields which acquire masses from SUSY breaking effects. Their masses are expected to be of order the gravitino mass but their couplings with the MSSM matter content are suppressed by a high energy scale . The moduli ($`\varphi `$) decay width will now be parameterized now as

$$\mathrm{\Gamma }_\varphi =\frac{1}{2\pi }\frac{m_\varphi ^3}{M_I^2}$$
(11)
where $`M_I`$ is the unification scale ($`10^{12}`$ GeV) which acts as an effective suppression scale.

The reheating temperature is obtained using eqs. (6) and (11)

$`{\displaystyle \frac{T_{RH}}{5MeV}}={\displaystyle \frac{10^{14}GeV}{M_I}}\left({\displaystyle \frac{m_\varphi }{100GeV}}\right)^{3/2}\left({\displaystyle \frac{10.75}{g_{}(T_{RH})}}\right)^{1/4},`$ (12)
where, for example, $`g_{}=10.75`$ for $`T`$ $`𝒪`$($`110`$) MeV. This reheating temperature is shown as a function of the modulus mass in Fig. 1 for different values of $`M_I`$. The request that the modulus mass is larger than $`100500`$ GeV in order to allow for kinematical decays into neutralinos of suitable mass $`m_{\stackrel{~}{\chi }_1^0}50200`$ GeV, limits in practice the reheating temperature to be above $`3`$ GeV for the lowest scale on consideration $`M_I=10^{11}`$ GeV. This has important consequences for the relic density computation. As discussed in the introduction of this section, we need a temperature smaller than the typical $`T_Fm_{\stackrel{~}{\chi }_1^0}/20310`$ GeV, in order to increase the relic neutralino density. This implies that the scale $`M_I10^{11}`$ GeV is in the border of validity. On the other hand, for larger values we can obtain very easily interesting reheating temperatures. For example, for $`M_I=10^{12}`$ GeV we have $`T_{RH}>0.3`$ GeV. For the highest scale with an interesting phenomenological value of the neutralino–nucleus cross section, in the case of universality and moderate $`\mathrm{tan}\beta `$ , $`M_I=10^{14}`$ GeV, the lowest value of the reheating temperature corresponds to $`T_{RH}6`$ MeV. Larger values of the scale, $`M_I>10^{15}`$ GeV, producing also a large cross section, are possible in D-brane scenarios since non-universality in soft terms is generically present . In this case the constraint $`T_{RH}>`$ 1 MeV from nucleosynthesis can be translated into a constraint on the modulus mass $`m_\varphi >140`$ GeV.

When considering the direct decay of the moduli fields into LSP (and therefore, an additional quantity of non thermal LSP production) the evolution of the cosmological LSP abundance becomes more complicated. We have to consider now the coupled Boltzmann equations for the LSP, the moduli field and the radiation . The detailed equations used by us appear explicitly in Ref.(). These equations depend on $`m_\chi `$, the LSP mass, $`\overline{N}_{LSP}`$, the averaged number of LSP particles produced in the decay of the modulus field.

Let us discuss qualitatively the solution following the arguments used in ref. . For a $`T_{RH}`$ higher than $`T_F`$ the relic density will roughly reproduce the usual result given by eq. (4). However, for the interesting case for us when $`T_{RH}`$ is lower than $`T_F`$, neutralinos produced from modulus decay are never in chemical equilibrium, unlike the thermal production case reviewed in Section 2. As a consequence, its number density always decreases through pair annihilation. When the annihilation rate $`\sigma _{\stackrel{~}{\chi }_1^0}^{ann}vn_{\stackrel{~}{\chi }_1^0}`$ drops below the expansion rate of the universe, $`H`$, the neutralino freezes out. Then the relic density can be estimated as

$`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2`$ $`=`$ $`{\displaystyle \frac{3m_{\stackrel{~}{\chi }_1^0}\mathrm{\Gamma }_\varphi }{2(2\pi ^2/45)g_{}T_{RH}^3\sigma _{\stackrel{~}{\chi }_1^0}^{ann}v}}{\displaystyle \frac{h^2}{\rho _c/s_0}}.`$ (13)
This result is valid when there is a large number of neutralinos produced by the modulus decay. When the number is insufficient, they do not annihilate and therefore all the neutralinos survive. The result in this case is given by

$`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2`$ $`=`$ $`{\displaystyle \frac{3\overline{N}_{\stackrel{~}{\chi }_1^0}m_\chi \mathrm{\Gamma }_\varphi ^2M_P^2}{(2\pi ^2/45)g_{}T_{RH}^3m_\varphi }}{\displaystyle \frac{h^2}{\rho _c/s_0}}.`$ (14)
The actual relic density is estimated  as the minimum of (13) and (14).

Now we can apply the above equations to our case with intermediate scales. Using eqs. (11) and (12), we can write expressions (13) and (14) as

$`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2{\displaystyle \frac{M_Im_{\stackrel{~}{\chi }_1^0}}{\sigma _{\stackrel{~}{\chi }_1^0}^{ann}vm_\varphi ^{3/2}g_{}^{1/4}}}`$ (15)
$`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2\overline{N}_{\stackrel{~}{\chi }_1^0}m_{\stackrel{~}{\chi }_1^0}{\displaystyle \frac{m_\varphi ^{1/2}}{g_{}^{1/4}M_I}}`$ (16)
From these equations we can see that even with a large annihilation cross section, $`\sigma _{\stackrel{~}{\chi }_1^0}^{ann}v10^8`$ GeV-2, we are able to obtain the cosmologically interesting value $`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^21`$. For example for $`M_I10^{13}`$ GeV we obtain it when $`\overline{N}_{\stackrel{~}{\chi }_1^0}1`$, using eq. (15).

Previous expressions are useful but approximate. We have solved numerically the Boltzmann equations in several cases. In Fig. 2 we show in detail the results of the numerical computation. We have used a particular expression for the annihilation cross section introduced which is explained below. In this figure, the contours of constant relic neutralino density $`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2`$ as a function of $`m_\varphi `$ and $`\overline{N}_{\stackrel{~}{\chi }_1^0}`$ are shown, for fixed values of $`M_I`$ and $`m_{\stackrel{~}{\chi }_1^0}`$. In particular, we consider the cases $`M_I=10^{12},10^{13},10^{14}`$ GeV, with $`m_{\stackrel{~}{\chi }_1^0}=100`$ GeV. The corresponding reheating temperatures can be obtained from Fig. 1. Note that whereas many values of $`\overline{N}_{\stackrel{~}{\chi }_1^0}`$ correspond to a satisfactory relic density for $`M_I=10^{12}10^{13}`$ GeV, for the case $`M_I=10^{14}`$ GeV only a small range works.

Let us finally finish this section with a few words about the neutralino annihilation cross section. As it is well known, in most of the parameter space of the MSSM the neutralino is mainly pure bino, and as a consequence it will mainly annihilate into lepton pairs through $`t`$–channel exchange of right–handed sleptons. The $`p`$–wave dominant cross section is given by

$`\sigma _{\stackrel{~}{\chi }_1^0}^{ann}v8\pi \alpha ^2{\displaystyle \frac{1}{m_{\stackrel{~}{\chi }_1^0}^2}}{\displaystyle \frac{1}{\left(1+x_{\stackrel{~}{l}_R}\right)^2}}v^2,`$ (17)
where $`x_{\stackrel{~}{l}_R}m_{\stackrel{~}{l}_R}^2/m_{\stackrel{~}{\chi }_1^0}^2`$ and $`\alpha ^{}`$ is the coupling constant for the $`U(1)_Y`$ interaction. Taking $`m_{\stackrel{~}{l}_R}m_{\stackrel{~}{\chi }_1^0}100`$ GeV, $`\sigma _{\stackrel{~}{\chi }_1^0}^{ann}v`$ in eq. (17) becomes of the order of $`10^9`$ GeV-2 or smaller. Using eq. (4) an interesting relic abundance, $`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2>0.1`$, is obtained.

However, in Refs., we obtained that for $`M_I10^{12}10^{14}`$ GeV the neutralino LSP acquires a sizeable Higgsino component, as we commented before this is essential for enhancing the LSP-nucleon cross section. For a Higgsino-like LSP the dominant annihilation channel is the decay into a W boson pair. An upper bound on the cross section can be obtained :

$`\sigma _{ann}v`$ $`<=`$ $`{\displaystyle \frac{\pi \alpha _2^2}{2}}{\displaystyle \frac{1}{m_\chi ^2}}{\displaystyle \frac{(1x_W)^{3/2}}{(2x_W)^2}},`$ (18)
where $`\alpha _2`$ is the $`SU(2)`$ constant and $`x_W=m_W^2/m_\chi ^2`$. It is this last expression for the cross section the one used in the numerical example presented in Fig.1. For a LSP with $`m_\chi 100`$ GeV, the formula gives a maximum $`2\times 10^8`$ GeV-2. In the case that the LSP have sizeable components of both gaugino and Higgsino one expects smaller cross sections.

## 5 Conclusions and Comments

In conclusion we have analyzed the LSP relic density in supersymmetric models with intermediate unification scale. we find that the reheating temperature is of order 1 GeV ($`M_I10^{12}`$) and below down to the nucleosynthesis threshold. We have shown that in generic low $`T_{RH}`$ scenarios without non-thermal direct production is expected that the LSP relic density is too low. As an alternative, we have presented a concrete scenario where the reheating temperature is associated to the decay of oscillating moduli fields appearing in string theories, we show that the LSP relic density considerably increases with respect to the standard radiation-dominated case by the effect of the direct non-thermal production by the modulus field. The LSP can become a good dark matter candidate ($`0.010.1>\mathrm{\Omega }h^2<0.31`$) for $`M_I10^{12}10^{14}`$ GeV and $`m_\varphi 110`$ TeV.

Let us finally remark that the numerical value of $`\overline{N}_{\stackrel{~}{\chi }_1^0}`$ is in general model dependent. This was discussed in the context of supergravity in ref. . In this particular case both values $`\overline{N}_{\stackrel{~}{\chi }_1^0}1`$ and $`\overline{N}_{\stackrel{~}{\chi }_1^0}10^310^4`$ are plausible, depending on the characteristics of the supergravity theory under consideration.