text:
# Maximizing the entanglement of two mixed qubits

## Abstract

Two-qubit states occupy a large and relatively unexplored Hilbert space. Such states can be succinctly characterized by their degree of entanglement and purity. In this letter we investigate entangled mixed states and present a class of states that have the maximum amount of entanglement for a given linear entropy.

With the recent rapid developments in quantum information there has been a renewed interest in multiparticle quantum mechanics and entanglement. The properties of states between the pure, maximally-entangled, and completely mixed (separable) limits are not completely known and have not been fully characterized. The physically allowed degree of entanglement and mixture is a timely issue, given that entangled qubits are a critical resource in many quantum information applications (such as quantum computationDiVincenzo95 ; Vedral98 , quantum communicationSchumacher96 , quantum cryptographyEkert 91 ; Naik00  and teleportationBennett93 ; Bouwmeester97 ), and that entangled mixed states could be advantageous for certain quantum information situationsCleve99 .

The simplest non-trivial multiparticle system that can be investigated both theoretically and experimentally consists of two-qubits. A two-qubit system displays many of the paradoxical features of quantum mechanics such as superposition and entanglement. Extreme cases are well known and clear enough: maximally-entangled two particle states have been produced in a range of physical systems Aspect 82 ; Kwiat99 ; Turchette98 ; Turchette95 , while two-qubits have been encoded in product (non-entangled) states Braunstein99  via liquid NMR Chuang98 . Recently, however White et al. have experimentally generated polarization-entangled photons in both non-maximally entangled statesWhite99 , and general states with variable degree of mixture and entanglement White00 .

In this letter we explore theoretically the domain between pure, highly entangled states, and highly mixed, weakly entangled states. We will partially characterisepartial  such two-qubit states by their purity and degree of entanglement Bennett96 . Specifically, we address the question: What is the form of maximally-entangled mixed states, that is, states with the maximum amount of entanglement for a given degree of purity? Ishizaka et al. Ishizaka00  have proposed two-qubit mixed states in which the degree of entanglement cannot be increased further by any unitary operations (the Werner stateWerner89  is one such example). A numerical exploration of the entanglement - purity plane is used to establish an upper bound for the maximum amount of entanglement possible for a given purity, and vice versa. We derive an analytical form for this class of maximally-entangled mixed states (MEMS) and show it to be optimal for the entanglement and purity measures considered.

Currently a variety of measures are known for quantifying the degree of entanglement in a bipartite system. These include the entanglement of distillationBennett96 , the relative entropy of entanglementVedral98 , but the canonical measure of entanglement is called the entanglement of formation (EOF) Bennett96  and for an arbitrary two-qubit system is given byCoffman99

$`E_F(\widehat{\rho })=h\left({\displaystyle \frac{1+\sqrt{1\tau }}{2}}\right),`$ (1)
where $`h(x)=x\mathrm{log}_2(x)(1x)\mathrm{log}_2(1x)`$ is Shannon’s entropy function and $`\tau `$, the “tangle” Coffman99  (“concurrence” squared) is given by

$`\tau =𝒞^2=\left[\mathrm{max}\{\lambda _1\lambda _2\lambda _3\lambda _4,0\}\right]^2,`$ (2)
Here the $`\lambda `$’s are the square roots of the eigenvalues, in decreasing order, of the matrix, $`\widehat{\rho }\stackrel{~}{\widehat{\rho }}=\widehat{\rho }\sigma _y^A\sigma _y^B\widehat{\rho }^{}\sigma _y^A\sigma _y^B`$, where $`\widehat{\rho }^{}`$ denotes the complex conjugation of $`\widehat{\rho }`$ in the computational basis $`\{|00,|01,|10,|11\}`$, and is an anti-unitary operation. Since the entanglement of formation $`E_F`$ is a strictly monotonic function of $`\tau `$, the maximum of $`\tau `$ corresponds to the maximum of $`E_F`$. Thus in this letter we use the tangle directly as our measure of entanglement. For a maximally-entangled pure state $`\tau =1`$, while for a unentangled state $`\tau =0`$.

There exist for the degree of mixture of a state a number of measures. These include the von Neumann entropy of a state, given by $`\mathrm{S}=\mathrm{Tr}[\widehat{\rho }\mathrm{ln}\widehat{\rho }]`$ vonNeumann55 , and the purity $`\mathrm{Tr}[\widehat{\rho }^2]`$. In this letter we use the linear entropy given byBose00

$`𝒮_L={\displaystyle \frac{4}{3}}\left\{1\mathrm{Tr}\left[\widehat{\rho }^2\right]\right\}.`$ (3)
which ranges from 0 (for a pure state) to 1 (for a maximally-mixed state). The linear entropy is generally a simpler quantity to calculate and hence its choice here.

Let us now examine our two-qubit states and the region they occupy in the tangle-linear entropy plane. We begin by randomly generating two million density matrices representing physical states, and determining their linear entropy and tangle. In Fig. (1a) we display a subset of these results for thirty thousand points. We see that quite a large region of this plane is filled with physically acceptable states (obviously a maximally-mixed, maximally-entangled, state is not possible). Zyczkowski et. al.Zyczkowski 99  have preformed similar numerical studies but their work focused on how many entangled states are in the set of all quantum states.

In Fig. (1a) we have also explicitly plotted the tangle versus linear entropy for the Werner state, a mixture of the maximally-entangled state and the maximally-mixed state Werner89 :

$`\widehat{\rho }={\displaystyle \frac{1\gamma }{4}}I_2I_2+\gamma |\mathrm{\Phi }_+\mathrm{\Phi }_+|,`$ (4)
where $`I_2`$ is the identity matrix and $`|\mathrm{\Phi }_+=\frac{1}{\sqrt{2}}\left[|0|0+|1|1\right]`$. We have labelled our orthogonal qubit states by $`|0`$ and $`|1`$. This Werner state is entangled (inseparable) for $`\gamma >1/3`$Bennett96a  and maximally-entangled when $`\gamma =1`$. The results from Fig. (1a) clearly indicate a class of states that have a larger degree of entanglement for a given linear entropy than the Werner states. We also generated a second set of data (by random perturbations about the maximally-entangled mixed states) so as to examine the boundary of possible states, which in the previous data set was a sparsely populated region. As can be seen in Fig. (1b), a definite boundary to the physically possible states exists.

Let us now analytically determine the form of these maximally-entangled mixed states (MEMS). As our starting point, let us consider the Werner state given by (4). How can one increase its degree of entanglement without changing its purity, or, alternatively, how can one increase its linear entropy given a certain degree of entanglement? In deriving our ansatz we will note the following points:

* In the Werner state (4) all the entanglement arises from the $`\gamma |\mathrm{\Phi }_+\mathrm{\Phi }_+|`$ term, and hence, to leave the degree of entanglement fixed while increasing the linear entropy this term needs to remain untouched. Local unitary operations will not affect the degree of entanglement or linear entropy.
* The $`I_2I_2`$ term of the Werner states comprises the maximally-mixed state. It can be written as an equal incoherent mixture of the four Bell states

$`|\mathrm{\Psi }_\pm `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|0|1\pm |1|0\right],`$ (5)
$`|\mathrm{\Phi }_\pm `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|0|0\pm |1|1\right],`$ (6)
In our ansatz, if we increase the amount of any of the $`|\mathrm{\Psi }_\pm `$ or $`|\mathrm{\Phi }_{}`$ Bell states, then the net entanglement in the total system generally decreases.
* In a general two-qubit density matrix there are two types of off-diagonal terms, those that represent the entanglement and those that represent single particle superposition. These single particle superposition terms can be set to zero by local linear operations, and so, by definition, cannot change the net entanglement.
* The diagonal elements of the two-qubit density matrix do not affect the system’s maximum entanglement (given a specified amount of $`|\mathrm{\Phi }_+\mathrm{\Phi }_+|`$. The diagonal elements, however, have a significant impact on the linear entropy.

These principles lead us to postulate an ansatz of the form

$`\widehat{\rho }=\left(\begin{array}{cccc}x+\frac{\gamma }{2}& 0& 0& \frac{\gamma }{2}\\ 0& a& 0& 0\\ 0& 0& b& 0\\ \frac{\gamma }{2}& 0& 0& y+\frac{\gamma }{2}\end{array}\right).`$ (11)
This comprises a mixture of the maximally-entangled Bell state $`|\mathrm{\Phi }_+`$ and a mixed diagonal state (who populations are specified by the real and non-negative parameters a,b,x,y). Without loss of generality we choose $`\gamma `$ to be a positive real number, which ensures that the ansatz density matrix is positive semi-definite. From normalization,

$`x+y+a+b+\gamma =1,`$ (12)
The linear entropy is simply given by

$`𝒮_L={\displaystyle \frac{4}{3}}\left\{1a^2b^2x^2y^2\gamma \left(x+y\right)\gamma ^2\right\}`$ (13)
with the concurrence given by

$`𝒞`$ $`=`$ $`\mathrm{Max}[\gamma 2\sqrt{ab},0].`$ (14)
To determine the form of the two-qubit maximally-entangled mixed states, we begin by specifying that the concurrence $`𝒞`$ must be greater than zero. Thus $`𝒞=\mathrm{Max}[\gamma 2\sqrt{ab},0]=\gamma 2\sqrt{ab}0`$ and therefore is maximized when $`𝒞=\gamma `$. This requires either $`a=0`$ and/or $`b=0`$ (without loss of generality we set $`b=0`$). Using the normalization constraint given by (12), the linear entropy is given by

$`𝒮_L={\displaystyle \frac{4}{3}}\left\{2a+\left(\gamma +2x\right)\left(1a\gamma \right)2x^22a^2\right\}.`$ (15)
Calculating the turning point of (15) we find that $`𝒮_L/x=0`$ when either $`x=0`$ (a minimum) or $`2x=1a\gamma `$ (a maximum) and $`𝒮_L/a=0`$ when either $`a=0`$ (a minimum) or $`4a=22x\gamma `$ (a maximum). First examining the $`𝒮_L/x`$ stationary solution and the maximum given by $`2x=1a\gamma `$, we observe that this condition requires $`x=y`$. If $`a=1\gamma `$ then the stationary point corresponds to a turning point. We now need to examine several parameter regimes to determine the optimal solution. The first region has concurrence values in the region $`2/3𝒞\gamma 1`$. In this region the optimal situation occurs when $`x=0`$ and $`a=1\gamma `$. This means the maximally-entangled mixed state has the form

$`\widehat{\rho }_{\mathrm{MEMS}}=\left(\begin{array}{cccc}\gamma /2& 0& 0& \gamma /2\\ 0& 1\gamma & 0& 0\\ 0& 0& 0& 0\\ \gamma /2& 0& 0& \gamma /2\end{array}\right).`$ (20)
The second regime occurs for $`0𝒞\gamma 2/3`$. In this case the optimal solution occurs when $`a=1/3`$ and $`x+\gamma /2=1/3`$. The optimal maximally-entangled mixed state in this region has the form

$`\widehat{\rho }_{\mathrm{MEMS}}=\left(\begin{array}{cccc}1/3& 0& 0& \gamma /2\\ 0& 1/3& 0& 0\\ 0& 0& 0& 0\\ \gamma /2& 0& 0& 1/3\end{array}\right).`$ (25)
In this case the diagonal elements do not vary with $`\gamma `$. Combining both these solutions, we can obtain (up to local unitary transformations) the following single explicit form for the maximal entangled mixed state:

$`\widehat{\rho }_{\mathrm{MEMS}}=\left(\begin{array}{cccc}g(\gamma )& 0& 0& \gamma /2\\ 0& 12g(\gamma )& 0& 0\\ 0& 0& 0& 0\\ \gamma /2& 0& 0& g(\gamma )\end{array}\right),`$ (30)
where

$`g(\gamma )=\{\begin{array}{cccc}\gamma /2& & & 𝒞\gamma 2/3\\ 1/3& & & 𝒞\gamma <2/3\end{array}.`$ (33)
The degree of entanglement for this maximally-entangled mixed state is simply $`\tau =\gamma ^2`$, while the linear entropy has the form

$`𝒮_L`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left[4g(\gamma )\left(23g(\gamma )\right)\gamma ^2\right].`$ (34)
In figure (1) we have plotted the tangle versus the linear entropy for the Werner state, and the numerically determined maximally-entangled mixed state. Our analytic expression for the state (30) perfectly overlays the numerically generated optimal curve. It is clear that these states have a significantly greater degree of entanglement for a given linear entropy than the corresponding Werner states. The maximally-entangled mixed state and Werner state curves join each other at two points in the tangle - linear entropy plane. The first and most obvious point occurs at $`(\tau ,S_L)=(1,0)`$ (here both states are maximally entangled). The second point occurs at $`(\tau ,S_L)=(0,8/9)`$. Here the two states are given by,

$`\widehat{\rho }_{\mathrm{WERNER}}=\left(\begin{array}{cccc}\frac{1}{3}& 0& 0& \frac{1}{6}\\ 0& \frac{1}{6}& 0& 0\\ 0& 0& \frac{1}{6}& 0\\ \frac{1}{6}& 0& 0& \frac{1}{3}\end{array}\right),\widehat{\rho }_{\mathrm{MEMS}}=\left(\begin{array}{cccc}\frac{1}{3}& 0& 0& 0\\ 0& \frac{1}{3}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{3}\end{array}\right).`$ (43)
Neither state is entangled. We observe that $`\widehat{\rho }_{\mathrm{MEMS}}`$ at this point has no nonzero off-diagonal elements but the Werner state does. The maximally-entangled mixed state is entangled as soon as the off-diagonal elements are nonzero ($`\gamma >0`$, while the Werner state requires $`\gamma >1/3`$ to be entangled). Though $`\widehat{\rho }_{\mathrm{WERNER}}`$ and $`\widehat{\rho }_{\mathrm{MEMS}}`$ have different forms they have the same degree of entanglement (zero) and linear entropy. Because of the way the maximally-entangled mixed state has been constructed, it never attains a linear entropy $`S_L=1`$. The Werner state attains this point because of its maximally-mixed component.

To confirm that our analytic solution is optimal and that no density matrix has a greater degree of entanglement for a given linear entropy than the state (30), we generated one million further random density matrices. We found that the maximally-entangled mixed state is indeed optimal. It is interesting to note, however, that the state is only optimal for mixture measures based on $`\mathrm{Tr}\left[\widehat{\rho }^2\right]`$; if instead the degree of mixture is measured for instance by the entropyvonNeumann55  the state is not optimal.

Lastly, how does our class of maximally-entangled mixed states compare with those predicted by Ishizaka and HiroshimaIshizaka00 ? Ishizaka’s two-qubit mixed states, the Werner state being a specific example, were chosen so that the degree of entanglement of such states cannot be increased further by unitary operations. In contrast we have derived a class of states that have the maximum amount of entanglement for a given linear entropy (and vice versa). Therefore our states are members of the Ishizaka et al. class by definition, although they were not explicitly consideredIshizaka00 . The Ishizaka et al. result indicates that a maximally-entangled mixed state cannot have the degree of entanglement increased by unitary operations. This state can however have its entanglement increased by a simple and experimentally realizable non-unitary concentration protocol recently proposed by Thew and MunroThew00 . Such a protocol is a based on generalization of the Procrustean method originally introduced for pure statesBennett96b  and recently demonstrated experimentallyKwiat 00 . In Figure (2) we display the results of the concentration protocol for two initial conditions. The solid curves represent a range of states that are obtainable, from the maximally-entangled mixed state, as the concentration protocol is applied to improve the output state characteristics. We observe that for all $`\gamma `$ the output characteristics can be significantly improved (solid grey lines). In fact for $`\gamma 2/3`$ the maximally-entangled mixed state can be concentrated up the dashed curve to a maximally-entangled pure state.

To summarize, we have discovered a class of partially entangled mixed two-qubit states that have the maximum amount of entanglement for a given linear entropy. An analytical form for these states was derived and they were shown to have significantly more entanglement for a given degree of purity than the Werner states. The properties of these states are still largely unknown and require significant exploration. Open questions such as “can such states be realized experimentally”, “to what extent do they violate Bell inequalities?”, and “do they have information processing advantages over other states” are the subject of current investigation.

We wish to thank K. Nemoto and G. J. Milburn for encouraging discussions. WJM and AGW would like to acknowledge the support of the Australian Research Council while DFVJ would like to thank the University of Queensland for their hospitality during his visit.


format: text

text:
# GEST

## Abstract

Galactic Exoplanet Survey Telescope (GEST) was proposed for a discovery mission to search for microlensing terrestrial planets toward the Galactic bulge and also Kuiper Belt Objects (KBOs) that are believed to hold vital information of the early history of the solar system. Microlensing planet search method is hinged on photometric singular behavior of lensing (refraction) that is due to discontinuities in the distribution of photon paths whose phenomenon is better known in total reflection. The singularities (caustics) directly translate into potentially large planetary signals but the small size of the caustics imposes four essential requirements for earth mass planet searches: massive survey, angular resolution, temporal resolution, and continuous monitoring capability in a statistically stable observing condition. A 2m scale wide FOV space telescope with a large focal plane such as the GEST meets the needs.

Habitability of earth mass planets in the habitable (liquid water) zone is suspected to depend crucially on a giant planet near the snowline. The large mass of the giant planets makes the planetary caustics larger and easier to detect them in abundance. In microlensing, one needs not wait 12 years to detect a jupiter at the Jupiter orbit (5 AU) unlike in other planet search methods (doppler, astrometry, transit) because the planetary microlensing events will complete their courses within 50 days or so – 70 times shorter!

Close-in giant planets of the bulge stars will be registered as transits in the microlensing data base and will form a valuable platform for comparison studies of the planets from different detection methods and different environments.

When the Galactic bulge is behind the sun, the GEST will be ideal for high resolution mapping of the large scale structures which will include weak lensing by dark matter, quasar lensing and host galaxies, strong lensing by clusters, and a large volume of galaxies, which with selective follow-up measurements of whose redshifts and IR images will shed light on the dark stuff (matter, energy, essence, extra dimensions, topological defects, … ).

Gravitation, Relativities, Microlensing and Planets Microlensing involves both special and general relativities in a most trivial way: the spacial warp due to lensing masses is scanned in time by specical relativistic particles (photons) emitted from the lensed source and the general relativistic effect is contained in a factor 2. In microlensing, none of the 3 D’s of lensing – delay, deflection, and distortion – are measured unlike in cosmological lensing because the time of flight differcences and the image sizes and separations are too small. What is measured is the total magnification of the images that varies in time because of the relative motion of the lens and the lensed as seen from the observer. Microlensing light curves are essentially smooth – characteristically smooth – except at caustic crossing discontinuities. Even the discontinuities are characteristically orderly and smoothed over the size of the source (lensed) star. So, given sufficient quality of data, microlensing light curves are easily distinguished and identified against the background of variable stars and flare stars.

A microlensing planetary system will be discovered as a low multiplicity $`n`$-point lens. The interference pattern of the Newtonian potentials of the masses of the host star(s) and planets determines the angular positions, sizes, and strengths of the caustics – the planetary signal generators. The caustic curve of a planetary system lens consists of a stellar caustic and planetary caustics. See figure 1 for the case of planetary binary lenses. The stellar caustic dominates the behavior of the planetary system as a gravitational lens except in the planetary caustic regions. The planetary caustics behave like small magnifying glasses that modify the would-be single lens light curve by the host star. The “size of the planetary caustic” (circles in figure 1) of an earth mass planet with mass fraction $`ϵ=3\times 10^6`$ is $``$ a few $`\mu `$as.

Ground-based Searches for Earth Mass Planets? The previous estimations of the feasibility of finding earth mass planets from the ground were based on the idea that turn-off stars (brighter than main sequence stars and smaller than giant stars) can be surveyed as source stars tytler ; emplanet ; gould ; sackett . The caustic regions in figure 1 were calculated with the size of a turn-off star ($`3R_{}`$). However, it is a false assumption that was derived from a color magnitude diagram of ground-based observations of the crowded Galactic bulge fields. The alleged turn-off stars are mostly blended main sequence stars (see figure 3). More recent HST data show that turn-off stars are indeed rare holtz  as one would expect from stellar evolutions. This is detrimental to ground-based searches for earth mass planets because of the large seeing disks. It is also the case that the moon closely accompanies the bulge during the high bulge season. With a 4m scale dedicated telesope at the best site (such as dedicated VISTA in Paranal), one can detect earth mass planets, but it takes a very long time to acquire sufficient statistics. The details are being calculated and will be reported in the revised version of gest .

GEST, Planets, KBOs, and Large Scale Structures: GEST is an ongoing effort to search from space for microlensing planets in the Galactic bulge and the disk. The lack of close-in giants in the globular cluster 47 Tucanae indicates the possibility of metallicity constraints on the formation of planets or density constraints on the survival of the planets gilliland . Galactic bulge is relatively metal rich ($``$ solar metallitcy in average) and is expected to be rich with terrestrial planets which may have been habitable when the Galaxy was younger. Gonzalez et al. GHZ  predicts massive terrestrial planets in the bugle, and GEST can explore the metallicity effect on the terrestrial planet mass limit. Ejection of planets is a commonality according to numerical studies, and GEST can detect the free-floating planets further constraining the planet formation and evolution scenarios. Microlensing is the only means to find mature free-floating planets.

Recent studies of habitability of planetary systems find a strong correlation between habitable terrestrial planets and jovian planets lunine , and this throws a double jeopardy for most of the planet search efforts: small mass of terrestrial planets and long orbital period of jovian planets. Microlensing is a peculiar method that can handle the both challenges with ease.

In cosmological lensing, information is extracted from the spatial intensity distribution of the photons in the images. In order to reconstruct the mass distributions and extract the cosmological parameters accurately, it is crucial to find systems that offer sufficient number of constraints which in turn requires wide and deep imaging with an angular resolution of a small fraction of an arcsecond. Excellent systems can be followed up by the NGST or ground-based adpative instruements for further information.

We thank the GEST/Discovery proposal team members, especially the PI D. Bennett, and many other supporters of large format space imaging.


text:
# Hartree-Fock-Bogoliubov calculation of charge radii of Sn, Ba, Yb, and Pb isotopes

## I INTRODUCTION

Odd-even staggering of isotope shifts is a common phenomena observed in many isotopic chains . The results are interpreted as due to changes in the charge radii of these isotopes, i.e., charge radii of odd-neutron nuclei are smaller than the average radii of their even neighbors .

Probable mechanisms so far presented for odd-even staggering of isotope shifts are; core polarization by valence neutrons , and mean field effects of a three-body (or a density-dependent) pairing force . Core polarization effects were first investigated by Talmi  as a possible mechanism for producing odd-even staggering of isotope shifts. He could show a remarkably good fitting to changes in the charge radii of Ca isotopes, though it is not certain whether the shell model calculation with a realistic force can reproduce absolute magnitude of isotope shifts. An alternative mechanism was suggested by Zawischa et al. . They have demonstrated the importance of a density dependent pairing as well as mean field effects of the pairing on the odd-even staggering of isotope shifts. They could obtain good fitting to changes in the charge radii of Ca, Ba, Sn, and Pb isotopes.

In this paper we investigate in detail the latter possibility employing the HFB theory with axially-symmetric deformation. We assume the Skyrme SLy4 parameterization and the delta-force pairing with linear density dependence . We solve HFB equations in the canonical basis. Odd nuclei are treated in the blocking approximation.

This paper is organized as follows. In Section 2, we review some theoretical formalism useful for our purposes. In Section 3, we give numerical results and discussions on Sn, Ba, Yb, and Pb isotopes. Summary and conclusions are given in Section 4.

## II FORMALISM

We follow the HFB scheme prescribed by Reinhard et al . We start with the Hamiltonian,

$$H=T+V^{Skyrme}+V^{pair},$$
(1)
$$V^{pair}(1,2)=V_0\frac{1𝝈_1𝝈_2}{4}\delta (𝒓_1𝒓_2)\left(1\frac{\rho (\frac{𝒓_1+𝒓_2}{2})}{\rho _c}\right).$$
(2)
Here $`V^{Skyrme}`$ is the Skyrme force and $`V^{pair}`$ is the density-dependent delta-force pairing . The pairing force involves two parameters, i.e., $`V_0`$ sets the strength of the force and $`\rho _c`$ determines its spatial dependence.

The BCS ansatz for the wave function of a pairing many-body system

$$|\mathrm{\Phi }_{BCS}=\underset{i>0}{}(u_i+v_ia_{i}^{}{}_{}{}^{}a_{\overline{i}}^{}{}_{}{}^{})|$$
(3)
requires that the single particle states are orthonormalized and that occupations add up to the total particle number $`\mathrm{\Sigma }v_\alpha ^2=N`$. This ansatz carries two key densities . One is the one-body density,

$$\rho _q(𝒓)=\underset{\sigma }{}\mathrm{\Phi }_{BCS}|a_{𝒓\sigma q}^{}a_{𝒓\sigma q}|\mathrm{\Phi }_{BCS}=\underset{i}{}v_i^2\underset{\sigma }{}\varphi _i(𝒓,\sigma ,q)\varphi _i^{}(𝒓,\sigma ,q),$$
(4)
and the other is the pairing density,

$$\stackrel{~}{\rho }_q(𝒓)=\underset{\sigma }{}(2\sigma )\mathrm{\Phi }_{BCS}|a_{𝒓\sigma q}a_{𝒓\sigma q}|\mathrm{\Phi }_{BCS}=\underset{i}{}u_iv_i\underset{\sigma }{}\varphi _i(𝒓,\sigma ,q)\varphi _i^{}(𝒓,\sigma ,q).$$
(5)
Here $`𝒓`$, $`\sigma `$, and q denote the space, spin, and charge state of a nucleon, respectively.

We assume the Skyrme force in the particle-hole channel, while the pairing force in the particle-particle channel. Thus the energy separates into a mean-field part and a pairing part as

$$E[\varphi ,v]=\mathrm{\Phi }_{BCS}|H|\mathrm{\Phi }_{BCS}=E_{HF}[\varphi ,v]+E_{pair}[\varphi ,v].$$
(6)
$`E_{HF}`$ is the Skyrme Hartree-Fock energy and $`E_{pair}`$ is the pairing energy,

$$E_{pair}[\varphi ,v]=\frac{V_0}{4}\underset{q}{}\stackrel{~}{\rho }_q^2(𝒓)\left(1\frac{\rho (𝒓)}{\rho _c}\right)𝑑𝒓.$$
(7)
Variation of the energy E with respect to a single-particle wave function is written as

$$\frac{\delta E}{\delta \varphi _i^{}}=\widehat{}_i\varphi _i,$$
(8)
$$\widehat{}_i=v_i^2(\widehat{h}_{Skyme}+\mathrm{\Gamma }_{pair}(𝒓))+u_iv_i\mathrm{\Delta }_q(𝒓).$$
(9)
Variation of $`E_{HF}`$ yields the Skyrme Hartree-Fock hamiltonian $`\widehat{h}_{Skyrme}`$ , while variation of $`E_{pair}`$ yields the pair mean-field $`\mathrm{\Gamma }_{pair}`$as well as the gap potential $`\mathrm{\Delta }_q`$. The pair mean-field $`\mathrm{\Gamma }_{pair}`$comes from the variation of the density $`\rho _q(𝒓)`$ in Eq. (7). It has the form

$$\mathrm{\Gamma }_{pair}(𝒓)=\frac{V_0}{4\rho _c}\left(\stackrel{~}{\rho }_{n}^{}{}_{}{}^{2}(𝒓)+\stackrel{~}{\rho }_{p}^{}{}_{}{}^{2}(𝒓)\right).$$
(10)
The gap potential $`\mathrm{\Delta }_q`$ comes from the variation of the pairing density $`\stackrel{~}{\rho }_q(𝒓)`$ in Eq. (7). It has the form

$$\mathrm{\Delta }_q(𝒓)=\frac{V_0}{2}\stackrel{~}{\rho }_q(𝒓)\left(1\frac{\rho (𝒓)}{\rho _c}\right).$$
(11)
It is noted that the pair mean-field $`\mathrm{\Gamma }_{pair}`$does not depend on the charge state. It is always repulsive and arises only from the density dependent pairing. On the other hand, the gap potential $`\mathrm{\Delta }_q`$ depends on the charge state q and arises from any pairing force irrespective of its density dependence. For a density parameter $`\rho _c`$ lower than the central density, the gap potential will be repulsive inside the nucleus and it is attractive at the nuclear surface.

A single particle hamiltonian $`_i`$ of Eq. (9) depends on the state i. It is not hermitian in contrast to HF + BCS calculations. Therefore variation actually needs to take care of the orthonormality of the wave functions as

$$\frac{\delta }{\delta \varphi _i^{}}\left(E\underset{ij}{}\mathrm{\Lambda }_{ij}\left(\varphi _j^{}(\xi )\varphi _i(\xi )𝑑\xi \delta _{ij}\right)\right)=0,$$
(12)
where $`\xi `$ represents all variables. The equation leads to a generalized mean-field equation,

$$\widehat{}_i\varphi _i=\underset{j}{}\mathrm{\Lambda }_{ji}\varphi _j.$$
(13)
Lagrange multipliers $`\mathrm{\Lambda }_{ij}`$ constitute a symmetric matrix $`\mathrm{\Lambda }_{ij}=\mathrm{\Lambda }_{ji}`$ .

Variation of the energy with respect to the occupation probability under the particle number constraint yields

$$\frac{d}{dv_i}\left(E\lambda \underset{j}{}v_j^2\right)=2v_ih_{ii}+(\mathrm{\Delta }_{ii})_q\left(\frac{v_i^2}{u_i}u_i\right)=0,$$
(14)
$$h_{ii}=\varphi _i|\widehat{h}_{Skyrme}+\mathrm{\Gamma }_{pair}(𝒓)\lambda |\varphi _i,$$
(15)
and the pairing gap is calculated as

$$(\mathrm{\Delta }_{ii})_q=\underset{\sigma }{}|\varphi _i(𝒓,\sigma ,q)|^2\mathrm{\Delta }_q(𝒓)d𝒓.$$
(16)
Eq. (14) can be solved in the standard manner and yields

$$\left\{\begin{array}{c}v_i\\ u_i\end{array}\right\}=\sqrt{\frac{1}{2}\frac{1}{2}\frac{h_{ii}\lambda }{\sqrt{(h_{ii}\lambda )^2+\mathrm{\Delta }_{ii}^{}{}_{}{}^{2}}}}.$$
(17)
Both Eqs. (13) and (14) constitute the HFB equation. Eqs. (13) and (14) are solved iteratively employing the gradient method,

$$|\varphi _i^{n+1}=𝒪\{|\varphi _i^n𝒟(\widehat{}_i^n|\varphi _i^n\underset{j}{}\mathrm{\Lambda }_{ji}|\varphi _j^n)\}.$$
(18)
The parameter $`𝒟`$ determines the size of the step and $`𝒪`$ means Gram-Schmidt orthonormalization. $`|\varphi _i^n`$ and $`\widehat{}_i^n`$ represent a single particle state and the mean field in the n-th iteration, respectively. Lagrange multipliers $`\mathrm{\Lambda }_{ij}`$ are determined by the orthonormality of the states. They are

$`\mathrm{\Lambda }_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\varphi _i^n|\widehat{}_i^n+\widehat{}_j^n|\varphi _j^n.`$ (19)
## III Results and discussions

We have investigated changes in the charge radii of Sn, Pb, Ba, and Yb isotopes because long isotopic chains are observed in these nuclei and because they lie in the spherical, transitional and well deformed regions of the nuclear chart, respectively. The parametrization SLy4 is employed for the Skyrme force. The parameters were specifically devised to reproduce correct incompressibility modulus for the symmetric nuclear matter and to improve the isospin property of the force away from the $`\beta `$-stability line.

These features of the SLy4 force are especially favorable to the calculation of charge radii along a long isotopic chain. However, the SLy4 paramaters were determined without the pairing correlation and therefore they would have to be modified when they are used in HFB calculations. With the pairing parameters discussed in the next paragraph, it turns out that pair mean-fields $`\mathrm{\Gamma }_{pair}`$contribute to HF + BCS results of Sn, Pb, Ba, and Yb isotopes by an order of ten keV and the gap potentials $`\mathrm{\Delta }`$ by an order of hundred keV. Gap potentials have non-negligible effects on the odd-even mass difference of isotopes. However, the effects may be included in the pairing by adjusting its parameters and therefore we assume the original SLy4 parametrization in the following analysis.

The depth and the density parameters of the pairing force are taken as $`V_0`$ = – 1250 MeV$``$fm3 and $`\rho _c`$ = 0.140 fm-3, respectively. They were obtained by fitting to one-neutron separation enegies as well as odd-even sttagering of charge radii of the above isotopes. We assumed the same set of pairing parameters for all isotopes. With a fixed value of $`\rho _c`$, we determined the depth $`V_0`$ by fitting to one-neutron separation energies. We repeated this process for several different values of $`\rho _c`$ and found the best value of $`\rho _c`$ by comparing odd-even staggerings of charge radii. With this procedure we obtained a lower density parameter $`\rho _c`$ = 0.140 fm-3 than the saturation density $`\rho _c`$ = 0.160 fm-3.

In practice, calculated charge-radii also depend on the active pairing space. In the present paper, we take into account all bound levels with energies lower than $`e_F`$ \+ 1$`\mathrm{}\omega `$, where $`e_F`$ is the Fermi energy and 1$`\mathrm{}\omega `$ = 41A-1/3 MeV. With this prescription, however, we have encountered a problem of convergence. In each step of iterations, some of the highest levels go in and out from the active pairing space, which give rise to an oscillation of the binding energy. Thus for the first few hundred iterations, we excluded the gap potential $`\mathrm{\Delta }_q`$ from the mean field hamiltonian and solved the equation by diagonalization. After having had a pretty good convergence, we fixed the number of single-particle levels of the active pairing space and solved the full mean-field equation by employing the gradient step.

We have also studied the cut-off prescription of Bonche et al. . With their prescription, however, we could not reproduce odd-even staggering in light Ba isotopes if we included the positive evergy levels in the active pairing space. Except for this difficulty, our cut-off prescription gives similar results to those of Bonche et al. .

We assumed $`𝒟=1/(E_{kin}+50)`$ MeV-1 as a size of the gradient step. The parameter $`𝒟`$ varies during the iterations and $`E_{kin}`$ is the kinetic energy of the previous iteration. We have also assumed that the iteration has converged if $`|\mathrm{\Delta }E|<`$ 10-4 MeV for the energy and $`|\mathrm{\Delta }r|<`$ 10-4 fm for the charge radius between two successive iterations. (We did not have a good solution for 125Sn, where the charge radius fulfilled $`|\mathrm{\Delta }r|<`$ 10-4 fm, but the binding energy oscillated with amplitude of $`|\mathrm{\Delta }E|`$ 2 $`\times `$ 10-3 MeV.)

Single particle states were expanded in an deformed oscillator basis with axially-symmetry . Eleven major shells were assumed in the calculation. Eleven major shells are not enough to produce correct binding energies of heavy nuclei. However, by comparing the present results with those of fifteen major shells, we have found that neutron separation energies and isotopic shifts were little affected by the size of the basis space. Odd neutron nuclei were treated in the blocking approximation . Odd neutron contributes to the density $`\rho `$, but it does not contribute to the pairing density $`\stackrel{~}{\rho }`$.

Fig. 1 shows calculated results of Sn isotopes along with experiment; (a) neutron separation energies, (b) isotopic shifts in the charge radii of Sn isotopes normalized to the nucleus 120Sn ($`\mathrm{\Delta }r_c^2=r_c^2(A)r_c^2(120)`$), and (c) isotopic shifts in the charge radii of Sn isotopes. Experiment is denoted by open circles, while calculation is denoted by crosses. As a guide for eyes, they are connected with solid lines and dashed lines, respectively. We also modified changes in the charge radii of Fig. 1(b) by subtracting an equivalent of the liquid-drop difference, $`r_{LD}^2(A)=3/5r_0^2A^{2/3}`$ and $`r_0`$ = 1.2 fm.

Figs. 1(a) and 1(b) show good agreement with experiment both for neutron separation energies and changes in the charge radii of Sn isotopes. Close examination of Fig. 1(b) reveals that charge radii of odd neutron nuclei are smaller than the average radii of their even neighbors. This is a well-known phenomena of odd-even staggering of isotope shifts. In Fig. 1(c) the phenomena is more clearly exhibited by taking the difference of charge radii between neighboring isotopes.

Odd-even staggering of isotope shifts has been a longstanding problem in nuclear physics. Only recently it is realized  that a density dependence seems to be necessary for the pairing and that mean field effects of the pairing seems to be responsible for this phonomena. Indeed, we cannot predict odd-even staggering of charge radii within a framework of ordinary HF + BCS theory, even if we employ the density dependent pairing in the calculation. We need to take into account mean fields of the pairing to explain the phenomena.

Fig. 2(a) shows the HF + BCS calculation without mean fields of the pairing. As is mentioned above, we could not reproduce odd-even staggering of isotope shifts even if we employed the density dependent pairing in the calculation. Fig. 2(b) shows the HF + BCS calculation with the pair mean-field $`\mathrm{\Gamma }_{pair}`$but without the neutron gap-potential $`\mathrm{\Delta }_n`$ . Odd-even staggering is now well reproduced except for heavy Sn isotopes. Fig. 2(c) shows the HFB calculation with the gap potential but without the pair mean-field. The figure shows that the gap potential is also responsible for the odd-even staggering. The potential works especially well for heavy Sn isotopes. Reminding that Fig. 1(c) is the full HFB calculation where the pair mean-field and the gap potential are both included in the hamiltonian, we may conclude that the pair mean-field as well as the gap potential are both responsible for the occurence of odd-even staggering of Sn isotopes.

Fig. 2(d) shows the pair mean-fields of 109Sn, 110Sn, and 111Sn. Fig. 2(e) shows the neutron gap potentials of the same nuclei. In these figures we show monopole radial shapes of the potentials in the Legendre expantion. From Fig. 2(d) we see that the pair mean-field is repulsive inside a nucleus. The repulsive potential causes a shallow HF field and thus the nucleus expands. Since an even nucleus has a stronger pair mean-field than its odd neighbors, the charge radius of the even nucleus expands more than its odd neighbors. The same is true for the gap potential. Fig. 2(e) shows that the gap potential is repulsive inside a nucleus while it is attractive at the nuclear surface. The potential too, expands the nucleus. The potential is also stronger for an even nucleus than its odd neighbors, leading to a larger charge radius for the even nucleus than its odd neighbors.

It is noted that the pair mean-field has a larger effect than the gap potential on the odd-even staggering of isotope shifts. This is because that the pair mean-field $`\mathrm{\Gamma }_{pair}`$exists in the proton HF field even for proton closed-shell nuclei. On the other hand, the proton gap-potential $`\mathrm{\Delta }_p`$ vanishes for the proton closed-shell nuclei. For such nuclei, the neutron gap-potential $`\mathrm{\Delta }_n`$ affects only indirectly the charge radii of isotopes via the HF fields.

Fig. 3(a) shows one-neutron separation energies of Pb isotopes. Fig. 3(b) shows isotopic shifts in the charge radii of Pb isotopes normalized to the nucleus 208Pb. Fig. 3(c) shows isotopic shifts in the charge radii of Pb isotopes. Good agreement with experiment was obtained both for neutron speration energies and isotopic shifts.

Pb isotopes are similar to Sn isotopes in that the proton shell is closed. However, in contrast to Sn isotopes, the isotopic chain of Pb extends beyond the neutron shell closure at N=126. At the shell closure we see in Fig. 3(b) an abrupt change of charge-radius differences, which is called a “kink”.

We have also investigated mean field effects of the pairing on the kink of Pb isotopes. Mean fields of the pairing get stronger as an isotope goes away from 208Pb and strong mean fields could give rise to a larger charge radius. If this is the case, changes in the charge radii of Pb isotopes would be small for nuclei lighter than 208Pb, while they would be large for nuclei heavier than 208Pb. Fig. 3(b) shows isotopic shifts of Pb nuclei from our calculation, the relativistic mean field (RMF) theory , and experiment. Our calculation could not predict large enough isotope shifts for nuclei beyond 208Pb. Indeed our calculation is very similar to the HF + BCS calculation without the mean field effects of the pairing. Mean fields of the pairing seem to have only marginal effects on the occurence of the kink around 208Pb.

The kink around 208Pb may be primarily due to the shell structures in this region of nuclei. The 2g9/2 neutron level lying close to the continuum seems to be the key to the occurence of the kink. Indeed the kink is well reproduced with the RMF calculation. The kink is also reproduced with a nonrelativistic Skyrme calculation by omitting the Fock term of the spin-orbit interaction or by assuming the isospin dependence in the spin-orbit energy functional.

We have so far studied mean field effects of the pairing on the charge radii of spherical nuclei. In the following we will examine mean field effects of the pairing on the charge radii of transitional and well deformed nuclei.

Fig. 4(a) shows neutron separation energies of Ba isotopes. Fig. 4(b) shows isotopic shifts in the charge radii of Ba isotopes. As in the case of Pb isotopes, the isotopic chain extends beyond the shell closure at N = 82 and exhibits a kink around 138Ba. Fig. 4(c) shows odd-even staggering of Ba isotopes. Agreement with experiment is in general good.

The kink in Ba isotopes may be due to changes in the nuclear deformation. Proton shells of Ba nuclei are not closed and therefore most Ba nuclei are soft against nuclear deformation. Fig. 4(d) shows calculated quadrupole moments which agree well with experiment . As seen from this figure, the deformation decreases as the neutron number increases toward 138Ba and then the deformation increases beyond 138Ba. Thus changes in the charge radii are small for nuclei lighter than 138Ba, while they are large for nuclei heavier than 138Ba. We see in Fig. 4(b) that our calculation reproduces the kink fairly well, though the deviation from experiment becomes appreciable in the region of light Ba isotopes.

At this point we wish to mention the reason why we have employed a small density parameter $`\rho _c`$ = 0.140 fm-3 for the pairing instead of the saturation density $`\rho _c`$ = 0.160 fm-3. Fig. 5(a) shows neutron separation energies of Ba isotopes calculated with $`\rho _c`$ = 0.160 fm-3 and $`V_0`$ = –1000 MeV$``$fm3. (A strength of $`V_0`$ = – 1000 MeV$``$fm3 is necessary in order to fit to neutron separation energies.) Fig. 5(b) shows isotopic shifts between neighboring isotopes calculated with the same set of parameters. By comparing Fig. 4(c) and Fig. 5(b), we see that with this set of parameters odd-even staggering of isotope shifts does not agree well with experiment. Although we have shown Ba isotopes as an example, we have obtained good results also for Sn and Pb isotopes with $`\rho _c`$ = 0.140 fm-3 and $`V_0`$ = –1250 MeV$``$fm3.

In Fig. 6, we show neutron separation energies of Yb isotopes. Agreement with experiment is good. Figs. 6(b) and 6(c) show changes in the charge radii of Yb isotopes. In this case, however, agreement with experiment is rather poor especially for nuclei lighter than 165Yb. Our calculation suggests that odd-even staggering in the charge radii of Yb isotopes might be related to the odd-even staggering of nuclear deformation of these nuclei. For nuclei lighter than 165Yb, a larger deformation is calculated for an even isotope than its odd neighbors. The difference is enough to explain the odd-even sttagering of charge radii of these nuclei. For nuclei heavier than 165Yb, however, a larger deformation is calculated for an odd nucleus than its even neighbors. Mean fields of the pairing and nuclear deformation have opposite effects on the charge radii, leading to a rather complicated behavier of isotopic shifts in this region of nuclei.

## IV Summary and Conclusion

We have investigated mean field effects of the pairing on the charge radii of Sn, Pb, Ba, and Yb isotopes. There are two mean fields arising from the pairing energy functional. One is the pair mean-field $`\mathrm{\Gamma }_{pair}`$which comes from the density dependent part of the pairing. The pair mean-field does not depend on the charge state and therefore the same form of the pair mean-field appears in the HF fields of protons and neutrons. In addition, the potential exists even for proton closed-shell nuclei. The other is the gap potential which arises from the variation of the pairing density in the pairing energy functional. Therefore the gap potential arises also from the density independent pairing. In contrast to the pair mean field, however, the proton gap potential vanishes for proton closed-shell nuclei.

For the pairing we have employed in the present calculations, the pair mean-field is repulsive inside a nucleus. The repulsive potential causes a shallow HF field and the nucleus expands. Since an even nucleus has a stronger pair mean-field than its odd neighbors, the charge radius of an even nucleus expands more than its odd neighbors. The same is true for the gap potential. It is repulsive inside a nucleus and attractive at the nuclear surface. The potential too, expands the nucleus. Like the pair mean field, the gap potential is stronger for an even nucleus than its odd neighbors. These facts seem to be the key ingredient in the occurence of odd-even staggering of isotopic shifts.

The present analyses have shown that the lower density parameter $`\rho _c`$=0.140 fm-3 than the saturation density $`\rho _c`$=0.160 fm-3 gives better agreement with experiment in fitting to the odd-even staggering of isotope shifts of Sn, Pb, Ba, and Yb. We have also shown that the kink observed in the charge radii of Ba isotopes may be due to changes in the nuclear deformation. The effects are enough to explain the kink around 138Ba.

## ACKNOWLEDMENTS

We would like to express our sincere appreciation to Ken-ichiro Arita for many fruitful and critical discussions. We would like to thank M. Yamagami and K. Matsuyanagi, and also members of Nagoya Nuclear-Structure Seminor for many illuminating discussions.


text:
# 1 Introduction

## 1 Introduction

In supergravity and (super)string models there exists a plethora of scalar and fermionic fields (we will loosely call them moduli $`\mathrm{\Phi }`$) with masses of the order of the weak scale and gravitational-strength couplings to ordinary matter. If produced in the early Universe, such quanta behave like nonrelativistic matter and decay at very late times, eventually dominating the energy of the Universe until it is too late for nucleosynthesis to occur. This poses a serious cosmological problem . Typical examples of these dangerous relics are the spin-3/2 gravitino (the supersymmetric partner of the graviton), the scalar moduli which parametrize supersymmetric flat directions in moduli space and seem almost ubiquitous in string theory, and their fermionic superpartners.

The slow decay rate of the $`\mathrm{\Phi }`$-particles is the essential source of cosmological problems because the decay products of these relics will destroy the 4He and D nuclei by photodissociation, and thus successful nucleosynthesis predictions . The most stringent bound comes from the resulting overproduction of D $`+`$ 3He; this requires that the $`\mathrm{\Phi }`$-abundance relative to the entropy density $`s`$ at the time of reheating after inflation  should satisfy

$$\frac{n_\mathrm{\Phi }}{s}\stackrel{<}{_{}}10^{12},$$
(1)
where the exact bound depends upon the mass $`m_\mathrm{\Phi }`$.

The dangerous moduli can be produced in the early Universe in a number of ways111Here we will study only production of moduli quanta and we do not consider the possibility of moduli oscillations around their minimum.. They may be generated by thermal scatterings in the plasma during the process of reheating after inflation when the energy density of the Universe gets converted into radiation. The number density $`n_\mathrm{\Phi }`$ at thermalization is estimated to be of the order of $`n_\mathrm{\Phi }/s10^2(T_r/M_{\mathrm{Pl}})`$, where $`T_r`$ is the final reheating temperature and $`M_{\mathrm{Pl}}=1.2\times 10^{19}`$ GeV is the Planck mass. This gives a stringent upper bound on the reheating temperature after inflation of about $`10^9`$ GeV .

The generation of the dangerous relics may also proceed through non-thermal effects. Scalar moduli fields may be copiously created by the classical gravitational effects on the vacuum state between the end of inflation and the beginning of the matter/radiation phase . In such a case, the upper bound on the reheating temperature is very dependent on the inflationary model, but it may be as low as 100 GeV. Another possibility is represented by the parametric excitation of dangerous relics occuring after inflation because of the rapid oscillations of the inflaton field(s). Gravitino production is an interesting example of this kind  since the non-thermal generation of gravitinos in the early Universe can be extremely efficient and overcome the thermal production by several orders of magnitude222 Recently, this statement has been criticized in ref. , where it was shown that, in a simple Polonyi model, gravitino production is suppressed. Undoubtedly, as emphasized in refs. , the non-thermal gravitino production depends sensitively on the assumptions on the inflationary model and its couplings to the supersymmetry-breaking sector. Non-thermal production of gravitinos is expected to be significant and to pose a cosmological problem if – after inflation – there is a strong mixing between the inflaton sector and the supersymmetry-breaking sector (of the present vacuum). As it will become clear in the present paper, this is often the case because, during the first few inflaton oscillations when the value of the inflaton field is close to $`M_{\mathrm{Pl}}`$, generic Planck mass-suppressed interactions can give a large contribution to such mixing..

The scope of this paper is to show that the abundant production of dangerous relics is an unavoidable consequence of a period of preheating  after inflation. The corresponding limits on the reheating temperature are often very tight and more severe than the bound of $`10^9`$ GeV coming from the production of moduli via thermal scatterings during reheating.

During the first stage of preheating nonperturbative quantum effects lead to an extremely effective dissipational dynamics and explosive particle production. Particles coupled to the inflaton field can be created in a broad parametric resonance with a fraction of the energy stored in the form of coherent inflaton oscillations at the end of inflation released after a few oscillations. This rapid transfer of energy density populates the Universe with quanta in highly nonthermal states which can be viewed either as classical waves travelling through the Universe  or as quantum particles in states with large occupation numbers. In the second stage of preheating, called semiclassical thermalization , rescatterings of the produced fluctuations smear out the resonance peaks in the power spectra and lead to a slowly evolving state in which the power spectra are smooth . The system begins to exibit a chaotic behaviour characteristic of a classical non-linear system with many degrees of freedom. In the course of the subsequent slow evolution, the power spectra propagate to larger momenta, eventually leading to a fully thermalized state.

In supersymmetric and superstring theories coupling constants and masses appearing in the Lagrangian have to be thought as functions of the dimensionless ratio $`\mathrm{\Phi }/M_{\mathrm{Pl}}`$. For instance, a generic coupling constant $`h`$ is in fact a function of the scalar moduli

$$h(\mathrm{\Phi })=h\left(\frac{\mathrm{\Phi }}{M_{\mathrm{Pl}}}\right)\left(1+c\frac{\delta \mathrm{\Phi }}{M_{\mathrm{Pl}}}+\mathrm{}\right),$$
(2)
where $`c`$ is a coefficient usually of order unity and $`\delta \mathrm{\Phi }=\mathrm{\Phi }\mathrm{\Phi }`$. This expansion introduces a gravitational coupling between the scalar moduli and the matter fields. If the underlying theory is supersymmetric, similar couplings of the scalar field $`\mathrm{\Phi }`$ and its fermionic superpartner to the matter fields will arise in the superpotential or in the Kähler potential.

Since moduli fields couple to any form of matter at least gravitationally, the production is present in a wide variety of inflationary models which have a preheating stage. In general, the equation of motion for the (quantum) field $`\delta \mathrm{\Phi }`$ during preheating is

$$\mathrm{}\delta \mathrm{\Phi }+m^2\delta \mathrm{\Phi }=\underset{i}{}h_iJ_i,$$
(3)
where the currents $`J_i`$ are functions of the (classical) fields created at preheating and have interaction vertices parametrized by the coupling constants $`h_i`$. In a cosmological setting particle creation in a time-varying classical background is usually caused by the time dependence of the effective mass of the field (consider, for instance, parametric resonance or gravitational creation of particles). As a result, particles are created in a squeezed state. On the contrary, creation of moduli is caused by the non-zero classical current and therefore they are created as coherent states. From a quantum mechanical point of view, this production mechanism is novel and differs from the more traditionally studied generation mechanism of squeezed states. Furthermore, the moduli production does not depend upon the particular properties of the modulus unlike the gravitino production  and it is true both for scalar and fermionic moduli.

In this paper we perform a calculation of moduli production during the preheating stage in various inflationary models. In sect. 2 we consider the cases of chaotic inflation with quartic and quadratic inflaton potentials, with and without couplings of the inflaton to a new scalar field $`X`$. The case of hybrid inflation is discussed in sect. 3, while sect. 4 contains a study of the production at preheating of modulinos, light fermions with only gravitational-strength couplings. In sect. 5 we summarize our results.

## 2 Moduli Production at Preheating in Chaotic Inflation

Let us first consider a simple model of chaotic inflation with quartic potential $`V(\varphi )=\frac{\lambda }{4}\varphi ^4`$. As discussed in sect. 1, we expect a generic coupling of the scalar moduli quanta $`\delta \mathrm{\Phi }`$ with the inflaton field $`\varphi `$ of the form $`\lambda J\delta \mathrm{\Phi }`$, where $`Jc\varphi ^4/(4M_{\mathrm{Pl}})`$. Let us denote the amplitude of the oscillating inflaton field as $`\varphi _0`$, while the frequency of the oscillations is $`\sqrt{\lambda }\varphi _0`$, where $`\varphi _0`$ is the zero-momentum mode of $`\varphi `$.

In this model fluctuations of the inflaton field do not grow during chaotization (they actually decrease due to the redshift) so the most efficient production of moduli fields take place at the end of the resonance phase. As we noted, one can describe the states we are considering either as a collection of interacting classical waves or as a collection of particles in modes with large occupation numbers. Therefore at late times it is possible to treat the system classically and to solve the relevant equations of motion on a lattice. The initial conditions for the classical problem are found solving the preceding quantum evolution by means of the appropriate Bogolyubov transformations (for more details see Ref. ).

The results of the numerical integration on a $`128^3`$ lattice for $`\lambda =10^{12}`$ are summarized in Fig. 1, where we plot the values of $`\varphi ^2`$, $`\mathrm{\Phi }^2`$ and of the inflaton zero mode squared $`\varphi _0^2`$ as a function of the conformal time  and in units $`M_{\mathrm{Pl}}^2`$. The field variance is given by $`\varphi ^2=[(2\pi )^3𝒱]^1d^3𝐤\varphi _𝐤^2`$, where $`𝒱`$ is the volume and $`\varphi _𝐤`$ is the Fourier $`𝐤`$-mode of the field $`\varphi `$. Here and in the following we have taken $`c=1`$. Time and particle momenta are given in units of $`\sqrt{\lambda }\varphi _0(0)`$, where $`\varphi _0(0)`$ is the initial value of the zero inflaton mode at the end of inflation. The time $`t=0`$, when inflation ends, is chosen as the moment of the first extremum of $`\varphi (t)`$. The maximum value of $`\varphi ^2`$ is achieved at the time $`t_r`$, which marks the end of the resonance phase, and, as shown in Fig. 1 is about $`10^7M_{\mathrm{Pl}}^2`$, corresponding to $`\mathrm{\Phi }^210^{15}M_{\mathrm{Pl}}^2`$. At the same time the zero mode is about $`\varphi _0(t_r)10^3M_{\mathrm{Pl}}`$.

We will now present some analytical estimates which approximately reproduce the numerical results of Fig. 1. To obtain the intensity of the moduli gravitational production which accompanies the creation and annihilation of the fluctuations generated during preheating, we first express the solution of the modulus equation of motion in terms of the retarded Green function

$$\mathrm{\Phi }(x)=\mathrm{\Phi }_{in}(x)+\lambda d^4x^{}G_{ret}(xx^{})J(x^{}).$$
(4)
Here $`\mathrm{\Phi }_{in}`$ is the initial configuration at time $`t=0`$ and the retarded Green function is given by

$$G_{ret}(xx^{})=\theta (tt^{})\frac{d^3𝐤}{(2\pi )^3}\frac{\mathrm{sin}\omega _k(tt^{})}{\omega _k}e^{i\stackrel{}{𝐤}\stackrel{}{𝐱}}.$$
(5)
In the case under consideration, the source term is given by $`J=c\varphi ^4/(4M_{\mathrm{Pl}})`$. Expanding the fields in Fourier modes, performing the integral over spatial coordinates, and taking $`\varphi _0\delta \varphi _𝐤`$ (justified by the numerical results of Fig. 1), we obtain

$$\mathrm{\Phi }_𝐤(t)=\frac{c\lambda }{M_{\mathrm{Pl}}}_0^t𝑑t^{}\frac{\mathrm{sin}\left[\omega _k(tt^{})\right]}{\omega _k}\varphi _0^3(t^{})\delta \varphi _𝐤(t^{}),$$
(6)
A crucial point is that the effect is linear in the inflaton fluctuations $`\delta \varphi _𝐤`$. Therefore, the phenomenon of moduli production is due to the oscillations (similar to oscillations of neutrino species in a medium) of the inflaton quanta into moduli in presence of the inflaton background $`\varphi _0`$, rather than to real rescatterings. From a pure quantum mechanical point of view, scalar moduli are generated as coherent states and not as squeezed states as in the usual parametric resonance.

The fluctuations $`\delta \varphi _𝐤`$ grow exponentially during the resonance phase as $`e^{\mu t}`$ where $`2\mu 0.07\sqrt{\lambda }\varphi _0(0)`$  and $`\varphi _0(0)0.3M_{\mathrm{Pl}}`$ is the initial value of the inflaton field at the end of inflation. Therefore the integral in eq. (6) can be approximately evaluated by assuming that the exponential determines the full time dependence of the integrand. Then we obtain that the maximum value of the variance of the modulus field, achieved at the end of the resonance phase, is given by

$$\mathrm{\Phi }^2\frac{c^2\lambda ^2\varphi _0^6(t_r)\varphi ^2}{\omega _k^2(t_r)\mu ^2M_{\mathrm{Pl}}^2}10^4c^2\frac{\varphi _0^4(t_r)}{M_{\mathrm{Pl}}^4}\varphi ^2.$$
(7)
Here we have used $`\omega _k(t_r)\sqrt{\lambda }\varphi _0(t_r)`$ as the typical frequency at the end of the resonance phase. Taking from Fig. 1 the values $`\varphi _010^3M_{\mathrm{Pl}}`$ and $`\varphi ^210^7M_{\mathrm{Pl}}^2`$ at $`t=t_r`$, eq. (7) gives $`\mathrm{\Phi }^2c^210^{15}M_{\mathrm{Pl}}^2`$, in very good agreement with the numerical result. Notice that the dependence on $`\lambda `$ dropped out from eq. (7). Our numerical calculation confirms that the maximum of $`\mathrm{\Phi }^2`$ is independent of the coupling $`\lambda `$ and it scales quadratically with $`c`$.

Since it is a good approximation to assume that the dominant particle production occurs around the time $`t_r`$, we find

$$n_\mathrm{\Phi }\omega _k(t_r)\mathrm{\Phi }^2\mathrm{and}n_\varphi \omega _k(t_r)\varphi ^2.$$
(8)
Using eq. (7) we obtain

$$\frac{n_\mathrm{\Phi }}{n_\varphi }10^8c^2\left[\frac{\varphi _0(t_r)}{10^3M_{\mathrm{Pl}}}\right]^4$$
(9)
at the end of the resonance phase.

The ratio of the number density of moduli in units of the entropy density $`s\rho ^{3/4}\lambda ^{3/4}\varphi _0^3(t_r)`$ is given by

$$\frac{n_\mathrm{\Phi }}{s}10^6c^2\left(\frac{10^{13}}{\lambda }\right)^{1/4}\left(\frac{\varphi _0^2(t_r)}{10^6M_{\mathrm{Pl}}^2}\right)\left(\frac{\varphi ^2}{10^7M_{\mathrm{Pl}}^2}\right).$$
(10)
For $`\lambda =10^{13}`$ (as fixed by the COBE normalization), this results contradicts the bound (1) by about six orders of magnitude (for $`c`$ of order unity) when the energy density in the scalar field is transferred to the energy density of a hot gas of relativistic particles. Notice that the ratio $`\frac{n_\mathrm{\Phi }}{s}`$ does not depend on the time of thermalization, because both $`n_\mathrm{\Phi }`$ and $`V^{3/4}`$ vary with the scale factor $`a`$ as $`a^3`$.

In this calculation (and throughout all this paper) we have made the implicit hypothesis that the moduli fields are light during the preheating stage. However, the same generic expansion in eq. (2) is expected to provide the moduli fields with a mass of the order of the Hubble rate $`H`$ during the preheating stage. Consider the simple model of chaotic inflation with potential $`V(\varphi )=\frac{\lambda }{4}\varphi ^4`$. The expansion of the coupling $`\lambda (\mathrm{\Phi }/M_{\mathrm{Pl}})`$ at the quadratic order in $`\mathrm{\Phi }/M_{\mathrm{Pl}}`$ leads to a term in the Lagrangian of the form $`\lambda (\delta \mathrm{\Phi }/M_{\mathrm{Pl}})^2\varphi ^4H^2(\delta \mathrm{\Phi })^2`$. Nevertheless, our numerical and analytical estimates are not affected by treating the moduli as light states. This is because moduli fields are generically produced with frequencies $`\omega `$ much larger than the Hubble rate ($`\omega \sqrt{\lambda }\varphi _0(t_r)H(t_r)`$ in the chaotic model $`\lambda \varphi ^4`$) and kinematically they may be regarded as massless. This does not come as a surprise since it is known that superheavy particles may be easily generated at preheating .

A simple variation of the model considered so far is the one in which a massless inflaton field is coupled to another field $`X`$, $`V=\frac{1}{4}\lambda \varphi ^4+\frac{1}{2}g^2\varphi ^2X^2`$. The oscillating field $`\varphi `$ amplifies fluctuations of the field $`X`$ via parametric resonance. The strength of the resonance can be parametrized by the $`q`$-parameter $`qg^2/4\lambda `$. Fluctuations of the inflaton field grow rapidly because of the back-reaction and the energy of the inflaton zero mode is quickly released into quanta of the $`X`$-field and the inflaton field.

For $`q=25`$, we numerically find that the maximum value of $`\varphi ^2`$ is about $`5\times 10^7M_{\mathrm{Pl}}^2`$ while $`X^2`$ is about one order of magnitude smaller, see Fig. 2. The resonance stage ends when $`\varphi _0^2(t_r)5\times 10^6M_{\mathrm{Pl}}^2\varphi ^2`$. If we again restrict ourselves to the interaction of the form $`c\lambda \delta \mathrm{\Phi }\varphi ^4/(4M_{\mathrm{Pl}})`$, moduli are produced through the scatterings of the inflaton quanta and – at the end of the resonance stage – $`\mathrm{\Phi }^210^{13}M_{\mathrm{Pl}}^2`$.

The analytic estimate in eq. (7) well reproduces this numerical result for $`\mathrm{\Phi }^2`$. The corresponding ratio $`\frac{n_\mathrm{\Phi }}{s}`$, derived from eq. (10), violates the bound (1) by seven orders of magnitude for $`c=1`$.

Another option is to have a coupling of the form $`g^2\delta \mathrm{\Phi }\varphi ^2X^2/(2M_{\mathrm{Pl}})`$. In such a case the production of moduli cannot take place through oscillations at the initial stage of the parametric resonance, but only via inverse decays of already created fluctuations of the $`X`$ field. This is equivalent to say that the effect in the Green function method is not linear in the fluctuations, but at most quadratic. However, moduli are produced at the initial stage of the semiclassical thermalization and the result is not much different from the previous case, see Fig. 3.

The main difference during the evolution arises because the effective couplings in the system depend upon the avereage value of the fluctuations of the moduli field. This changes effectively the parameter $`q`$ and the time development of the resonance which is very sensitive to any change of $`q`$.

We can now use the Green-function method to give an estimate which approximately reproduces the results found in Fig. 3. The Fourier modes of the moduli fields at time $`t`$ can be written as

$$\mathrm{\Phi }_𝐤(t)=\frac{cg^2}{2M_{\mathrm{Pl}}}_0^t𝑑t^{}\frac{\mathrm{sin}\left[\omega _k(tt^{})\right]}{\omega _k}\varphi _0^2(t^{})d^3pX_𝐩(t^{})X_{𝐤𝐩}(t^{}).$$
(11)
As expected, the moduli Fourier modes are quadratic in the fluctuations of the $`X`$ field. We can perform the time integral by assuming that the leading time dependence comes from the exponential growth of $`X_𝐩`$, which is proportional to $`e^{\mu _Xt}`$, with $`\mu _X0.1\sqrt{\lambda }\varphi _0(0)`$. Since the momentum distribution of the $`X`$ fluctuations is sharply peaked, the integral over $`𝐤`$ of the square of eq. (11) at the end of the resonance gives

$$\mathrm{\Phi }^2\left[\frac{g^2\varphi _0^2(t_r)X^2}{4M_{\mathrm{Pl}}\mu _X\omega _k(t_r)}\right]^210^3q^2\frac{\varphi _0^2(t_r)X^2^2}{M_{\mathrm{Pl}}^4}.$$
(12)
Here we have used $`\omega _k\sqrt{\lambda }\varphi _0`$ as the typical frequency. From Fig. 3 we gather that, at the maximum, $`X^210^7M_{\mathrm{Pl}}^2`$, corrisponding to $`\varphi _0^2(t_r)10^5M_{\mathrm{Pl}}^2`$. For $`q=25`$, eq. (12) gives $`\mathrm{\Phi }^210^{13}M_{\mathrm{Pl}}^2`$, in fair agreement with our numerical results.

In all cases we also have calculated the particle occupation numbers to make sure that the lattice spacing and the integration box are suitably chosen and the relevant particle momenta are well represented. Particle number densities in the comoving reference frame are shown in Fig. 4 for the same model which is represented in Fig. 3 and at the time $`\tau =180`$. Let us also note that the production of moduli during the preheating stage should be analyzed by letting the moduli be coupled to the matter fields in all the interaction terms. Therefore, our estimates provide a lower limit on the abundance of moduli states.

Reassured by the capability of reproducing the numerical results by analytical means, we have also considered the quadratic inflaton model with quartic coupling to another field $`X`$, $`V=\frac{1}{2}m^2\varphi ^2+\frac{1}{2}g^2\varphi ^2X^2`$. The $`q`$-parameter is $`q=g^2\varphi _0^2(0)/4m^2`$. It is also useful to introduce the redshifted resonance parameter at the end of the resonance stage $`q_r=q(\varphi _0(t_r)/\varphi _0(0))^2`$. Parametric resonance can fully develop in an expanding Universe if $`q_r\stackrel{>}{_{}}1`$ . Resonant production is most effective for fluctuations of fields which couple not too weakly but also not too strongly, those with $`q_r1`$. For $`q_r1`$, the maximal size of $`X`$ fluctuations is suppressed by non-linear effects . The zero mode rapidly drops at the end of the chaotization stage and all its energy is transferred to fluctuations ; the variance of the inflaton field is much larger than the value of the zero mode at the chaotization phase.

We now want to perform an analytical estimate of the number density of moduli fields for moderate $`q_r`$. If the mass term has a dependence on the modulus field $`\mathrm{\Phi }`$, $`m^2=m^2(\mathrm{\Phi }/M_{\mathrm{Pl}})`$, this will introduce a trilinear coupling $`m^2\varphi ^2\delta \mathrm{\Phi }/(2M_{\mathrm{Pl}})`$ between the quanta of the modulus and the inflaton quanta. In such a case, the Green-function method previously illustrated gives a modulus variance at the end of resonance equal to

$$\mathrm{\Phi }^2\frac{m^4\varphi ^2^2}{16\omega ^2(t_r)\mu ^2M_{\mathrm{Pl}}^2}\frac{\varphi ^2^2}{M_{\mathrm{Pl}}^2}.$$
(13)
Here we have taken $`\omega m`$ as the typical energy and $`\mu 0.2m`$ . For $`q_r1`$, which corresponds to $`q10^4`$ , the inflaton completely decays into fluctuations during the chaotization phase; at the end of it $`\varphi _0^2\varphi ^2n_\varphi /m10^6M_{\mathrm{Pl}}^2`$. Therefore, from eq. (13), we obtain $`n_\mathrm{\Phi }/n_\varphi 10^6`$.

Notice that we could obtain the same estimate using the Boltzmann equation for the number density $`n_\mathrm{\Phi }`$:

$$\dot{n}_\mathrm{\Phi }+3Hn_\mathrm{\Phi }\frac{m^4}{\omega ^4}\frac{n_\varphi ^2}{M_{\mathrm{Pl}}^2}.$$
(14)
Assuming again $`\omega m`$ and knowing that $`n_\varphi `$ grows exponentially with rate $`2\mu `$, we get $`n_\mathrm{\Phi }/n_\varphi n_\varphi /(\mu M_{\mathrm{Pl}}^2)`$.

At the end of the reheating stage, the energy density stored in the inflaton is converted into a relativistic thermal bath with temperature $`T_r(mn_\varphi )^{1/4}`$. The number density of moduli per entropy density will then be

$$\frac{n_\mathrm{\Phi }}{s}\frac{n_\mathrm{\Phi }}{n_\varphi }\frac{T_r}{m}10^6\frac{T_r}{m}.$$
(15)
Since the COBE normalization gives $`m10^{13}`$ GeV, the upper bound on the reheating temperature becomes

$$T_r\stackrel{<}{_{}}10^7\mathrm{GeV}.$$
(16)
If the coupling of the modulus is $`g^2\delta \mathrm{\Phi }\varphi ^2X^2/(2M_{\mathrm{Pl}})`$, the Green function method gives

$$\mathrm{\Phi }^2\left[\frac{g^2\varphi _0^2(t_r)X^2}{4M_{\mathrm{Pl}}\mu _X\omega (t_r)}\right]^2\left[\frac{q_rX^2}{0.2M_{\mathrm{Pl}}}\right]^2.$$
(17)
Here we have assumed $`\omega m`$ and $`\mu _X0.2m`$. Repeating the same argument used above, we find $`n_\mathrm{\Phi }/s10^7q_r^2X^2^2T_r/(M_{\mathrm{Pl}}^4m)`$. Since at the end of the resonance phase $`q_r1`$, we find that the bound (1) gives a limit on the reheating temperature of $`10^5`$ GeV if – for instance – we use $`q=10^6`$ and the maximum value of $`X^210^6M_{\mathrm{Pl}}^2`$ .

So far we have always considered moduli couplings arising from the scalar potential, but similar results are also obtained when the moduli couplings are derived from non-minimal kinetic terms. For instance, from an interaction of the form $`\delta \mathrm{\Phi }\dot{\varphi }^2/M_{\mathrm{Pl}}`$ we estimate

$$\mathrm{\Phi }^2\frac{\dot{\varphi }_0^2\dot{\varphi }^2}{M_{\mathrm{Pl}}^2\omega ^4}\frac{\varphi _0^2\varphi ^2}{M_{\mathrm{Pl}}^2}.$$
(18)
In the case of chaotic inflation with quadratic potential, eq. (18) is equivalent to eq. (13), valid for moduli coupled to the inflaton mass term.

## 3 Moduli Production at Preheating in Hybrid Inflation

In this section we study the moduli production during the preheating phase following a period of hybrid inflation . This scenario involves two scalar fields, the inflaton field $`\varphi `$, and the symmetry breaking field $`\sigma `$, and different mass scales and couplings. During inflation, the inflaton field $`\varphi `$ rolls down along a flat potential while the Higgs field $`\sigma `$ is stuck at the origin providing the vacuum energy density driving inflation. However, when $`\varphi `$ gets smaller than a critical value $`\varphi _c`$, both fields roll down very quickly towards their present minima, completing the inflationary phase.

We take the simplest hybrid inflation potential as suggested by Linde  333For other hybrid inflation models, including those motivated by supersymmetry, see ref. .

$$V(\varphi ,\sigma )=\frac{\lambda }{4}\left(v^2\sigma ^2\right)^2+\frac{1}{2}m_\varphi ^2\varphi ^2+\frac{1}{2}g^2\varphi ^2\sigma ^2.$$
(19)
This potential has a valley of minima at $`\sigma =0`$ for $`\varphi ^2>\varphi _c^2\frac{\lambda }{g^2}v^2`$. Most of the inflation occurs while $`\varphi `$ is slowly rolling down from its initial value to $`\varphi _c`$. During this phase $`\sigma =0`$ and the potential reduces to $`V=\frac{\lambda }{4}v^4+\frac{1}{2}m_\varphi ^2\varphi ^2`$. For $`\varphi <\varphi _c`$, the concavity of the potential in the $`\sigma `$ direction changes so that the $`\varphi `$-axis is a ridge of unstable maxima.

Inflation ends because of the growth of quantum fluctuations in $`\sigma `$. For $`\varphi >\varphi _c`$, these are quickly damped out, but beneath $`\varphi _c`$ they cause the fields to fall very quickly to global minima at $`\varphi =0`$ and $`\sigma ^2=v^2`$. There, they will execute damped oscillations.

The results of numerical integration for the case $`v=10^3M_{\mathrm{Pl}}`$ are shown in Fig. 5, where we have supposed that the modulus is coupled to the inflaton field and the Higgs field via the coupling $`g^2\mathrm{\Phi }\varphi ^2\sigma ^2/(2M_{\mathrm{Pl}})`$.

In hybrid models the coupling constant $`\lambda `$ does not have to obey the same restrictions as in the case of chaotic inflation. However, for the numerical integration we have chosen a rather low value, $`\lambda =10^8`$. This is because the particle number in the inflaton or Higgs excitations at the beginning of chaotization is $`n(k)1/\lambda `$ and to be able to treat moduli fields classically, we need occupation numbers $`n_\mathrm{\Phi }1`$. Since $`n_\mathrm{\Phi }n(k)`$, we need a rather small value of $`\lambda `$. For the run presented in Fig. 5 we have chosen $`g^2/\lambda =10`$.

The initial conditions were chosen in such a way that Higgs oscillations occur only around one of the minima, and domain creation and other effects leading to complete chaotisation during the first oscillation (see Ref. ) can be disregarded.

Our numerical findings show that the variance of the moduli fields does depend only upon the ratio $`g^2/\lambda `$.

The Green function method tells us that

$$\mathrm{\Phi }^2\frac{g^4}{M_{\mathrm{Pl}}^2}\frac{\varphi ^2^2}{\omega _\varphi ^2}\frac{\sigma ^2^2}{\omega _\sigma ^2}\frac{g^2}{\lambda }\frac{\varphi ^2^2}{M_{\mathrm{Pl}}^2},$$
(20)
since the typical frequencies are given by $`\omega _\varphi ^2\lambda \sigma ^2`$ and $`\omega _\sigma ^2g^2\sigma ^2`$. Using the same values of the couplings and mass scales adopted in Fig. 5 and reading from there the variance of the inflaton field $`\varphi ^210^2v^2`$, we get $`\mathrm{\Phi }^210^9v^2`$, in fair agreement with our numerical results. Notice that the numerical values of the inflaton and Higgs variances at their maximum are well reproduced by $`\varphi ^2(\lambda /g^2)v^2`$ and $`\sigma ^2v^2`$. Using these values we estimate $`\mathrm{\Phi }^2(\lambda /g^2)(v^4/M_{\mathrm{Pl}}^2)`$, in good agreement with our numerical results and with the fact that these results depend only upon the ratio $`g^2/\lambda `$.

If the final reheating takes place via decay of the lightest particles in the system, i.e. the Higgs particles with mass squared $`\lambda v^2`$, the final moduli number density-to-entropy-ratio is given by

$$\frac{n_\mathrm{\Phi }}{s}\frac{n_\mathrm{\Phi }}{n_\sigma }\frac{T_r}{\sqrt{\lambda }v}$$
(21)
which implies

$$T_r\stackrel{<}{_{}}\sqrt{\lambda }\left(\frac{10^3M_{\mathrm{Pl}}}{v}\right)\left(\frac{g^2/\lambda }{10}\right)\mathrm{\hspace{0.33em}10}^{11}\mathrm{GeV}.$$
(22)
This result illustrates how, in certain ranges of parameters, the coherent production of moduli can become more dangerous than the generation of moduli via thermal scatterings.

## 4 Modulino Production at Preheating

So far we have been considering the generation of scalar moduli during the preheating stage. In this section we consider the production of the fermionic partners $`\psi _\mathrm{\Phi }`$ of the moduli fields, sometimes dubbed modulini. As for the scalar moduli, the modulini are generically coupled to the matter fields through couplings suppressed by $`M_{\mathrm{Pl}}`$.

Consider – for instance – the superpotential $`W=\frac{\sqrt{\lambda }}{3}\varphi ^3`$ giving rise to the inflaton potential $`V=\lambda \varphi ^4`$. If the coupling $`\lambda `$ is a function of $`\mathrm{\Phi }/M_{\mathrm{Pl}}`$, where now $`\mathrm{\Phi }`$ has to be considered as a superfield, the Lagrangian will contain the following coupling

$$\frac{\sqrt{\lambda }}{M_{\mathrm{Pl}}}\overline{\psi }_\mathrm{\Phi }\psi _\varphi \varphi ^2,$$
(23)
where have denoted with $`\psi _\varphi `$ the inflatino, the fermionic partner of the inflaton field.

We can again make use of the Green function method to obtain the Fourier component of $`\psi _\mathrm{\Phi }`$. Using the retarded Green function for a fermionic field, we find

$$\psi _\mathrm{\Phi }(𝐤,t)=\frac{\sqrt{\lambda }}{M_{\mathrm{Pl}}}_0^t𝑑t^{}\left\{i\gamma ^0\mathrm{sin}\left[\omega _k(tt^{})\right]\frac{\stackrel{}{\gamma }\stackrel{}{k}}{\omega _k}\mathrm{cos}\left[\omega _k(tt^{})\right]\right\}\varphi _0^2(t^{})\psi _\varphi (𝐤,t^{}).$$
(24)
We can then approximatly estimate the quantity $`\overline{\psi }_\mathrm{\Phi }\psi _\mathrm{\Phi }`$ as

$$\overline{\psi }_\mathrm{\Phi }\psi _\mathrm{\Phi }\frac{\lambda }{M_{\mathrm{Pl}}^2}\frac{\varphi _0^4}{\omega ^2}\overline{\psi }_\varphi \psi _\varphi .$$
(25)
One can now use as a guide the recent results obtained in the theory of generation of Dirac fermions during and after inflation . During the first inflaton oscillations, the Fermi distribution function of the inflatini is rapidly saturated up to some maximum value of the momentum $`k`$, i.e $`n_k1`$ for $`k\stackrel{<}{_{}}k_{\mathrm{max}}`$ and it is zero otherwise. The resulting number density is therefore $`n_kk_{\mathrm{max}}^3`$. The value of $`k_{\mathrm{max}}`$ is expected to be roughly of the order of the inverse of the time-scale needed for the change of the mass of the inflatini. In the model under consideration, such a mass changes by an amount $`\sqrt{\lambda }\varphi _0`$ in a time scale $`(\sqrt{\lambda }\varphi _0)^1`$ and one expects $`k_{\mathrm{max}}\sqrt{\lambda }\varphi _0`$. Therefore

$$\overline{\psi }_\varphi \psi _\varphi \lambda ^{3/2}\varphi _0^3(0),$$
(26)
where we have used the fact that most of the fermion production takes place during the first oscillation. Under these circumstances, one can also assume that the number density of inflatini $`n_{\psi _\varphi }`$ and modulini $`n_{\psi _\mathrm{\Phi }}`$ is well measured by $`\overline{\psi }_\varphi \psi _\varphi `$ and $`\overline{\psi }_\mathrm{\Phi }\psi _\mathrm{\Phi }`$, respectively. This gives

$$\frac{n_{\psi _\mathrm{\Phi }}}{n_{\psi _\varphi }}\frac{\varphi _0^2(0)}{M_{\mathrm{Pl}}^2},$$
(27)
where we have taken $`\omega k_{\mathrm{max}}`$ in Eq. (25).

Finally, the ratio of the number density of modulini in units of the entropy density $`s\rho ^{3/4}\lambda ^{3/4}\varphi _0^3`$ is given by

$$\frac{n_{\psi _\mathrm{\Phi }}}{s}\left(\frac{\lambda }{10^{13}}\right)^{3/4}\left(\frac{\varphi _0(0)}{0.3M_{\mathrm{Pl}}^2}\right)^210^{11}.$$
(28)
This result will contradict the bound (1) by about an order of magnitude when the energy density in the scalar field is transferred to the energy density of a hot gas of relativistic particles.

Similar considerations apply to the inflaton model with quadratic potential. Suppose there is a coupling of the form $`m\varphi ^2\mathrm{\Phi }/M_{\mathrm{Pl}}`$ in the superpotential. This gives rise in the Lagrangian to a coupling of the modulini with the fermionic partners of the inflaton of the form $`\frac{m}{M_{\mathrm{Pl}}}\overline{\psi }_\mathrm{\Phi }\psi _\varphi \varphi `$. Using the Green function method we again obtain $`n_{\psi _\mathrm{\Phi }}n_{\psi _\varphi }(\varphi _0(0)/M_{\mathrm{Pl}})^2`$. This gives $`n_{\psi _\mathrm{\Phi }}/n_{\psi _\varphi }(m/M_{\mathrm{Pl}})^210^{12}`$ and the generation of dangerous relics is not large enough to provide a strong constraint on the reheating temperature.

In the case of the hybrid model though, the bound is more restrictive. In supersymmetric hybrid models, the coupling constants $`\lambda `$ and $`g`$ are related, and typically one finds $`\lambda /g^21`$ . The COBE normalization implies that typically $`v5\times 10^{15}`$ GeV. During the oscillations of the Higgs and the inflaton field, one expects that the number density of “Higgsini ” $`\psi _\sigma `$ and inflatini $`\psi _\varphi `$ is of the order of the number density of Higgs particles excited during the preheating phase, $`n_{\psi _\sigma }n_{\psi _\varphi }n_\sigma \lambda ^{3/2}v^3`$. If the modulus couples to the Higgs superfield via the interaction $`\sqrt{\lambda }(\mathrm{\Phi }/M_{\mathrm{Pl}})\sigma ^2\varphi `$ in the superpotential, the number density of modulini is expected to be

$$\frac{n_{\psi _\mathrm{\Phi }}}{n_\sigma }\left(\frac{v}{M_{\mathrm{Pl}}}\right)^2.$$
(29)
The corresponding bound on the reheating temperature can be obtained using the same arguments which led us to eq. (22); we find $`T_r\stackrel{<}{_{}}\sqrt{\lambda }\left(5\times 10^{15}\mathrm{GeV}/v\right)\mathrm{\hspace{0.33em}10}^{10}\mathrm{GeV}`$.

## 5 Conclusions

In this paper we have investigated the production of scalar and fermionic moduli fields which occurs during the preheating stage after inflation. If produced with too large abundances, the late decays of these relics may jeopardize the successful predictions of standard big-bang nucleosynthesis. Since the moduli couple to any form of matter at least gravitationally, their generation is unescapable in the large class of inflationary models which have a preheating stage. They are generated as coherent states initially. From this point of view, the production mechanism of moduli differs from the generation mechanism of squeezed states during the parametric resonance which has been studied in recent years. We have shown that moduli are efficiently produced during the preheating stage. The corresponding upper bound on the reheating temperature is model dependent and sensitive to the details of the stage of thermalization following the period of parametric resonance. It is often tighter than the bound of $`10^9`$ GeV found considering the generation of moduli via thermal scatterings during reheating.