text:
# Self-consistent fragmented excited states of trapped condensates.

## Abstract

Self-consistent excited states of condensates are solutions of the Gross-Pitaevskii (GP) equation and have been amply discussed in the literature and related to experiments. By introducing a more general mean-field which includes the GP one as a special case, we find a new class of self-consistent excited states. In these states macroscopic numbers of bosons reside in different one-particle functions, i.e., the states are fragmented. Still, a single chemical potential is associated with the condensate. A numerical example is presented, illustrating that the energies of the new, fragmented, states are much lower than those of the GP excited states, and that they are stable to variations of the particle number and shape of the trap potential.

Numerous properties of dilute Bose-Einstein condensates are well discussed by the Gross-Pitaevskii (GP) equation 1 ; 2 . For reviews see, e.g. 3 ; 4 ; 5 . If the condensate is trapped, this equation, which is equivalent to the so called nonlinear Schrödinger equation, possesses a discrete spectrum of stationary states. Being a nonlinear equation which has to be solved self-consistently, the solutions of the GP equation are generally called self-consistent states. The properties of GP self-consistent excited states and their formation have been amply discussed, see, e.g. 6 ; 7 ; 8 ; 9 ; 10 ; 11 ; 12 ; 13 ; 14 . The GP equation supports many well-known self-consistent excitations, such as vortex states 15 ; 16 ; 17 ; 18 ; 19 , and bright and dark solitons 20 ; 21 ; 22 ; 23 ; 24 . The observation of solitons in trapped Bose gases 25 ; 26  provides a striking manifestation of nonlinear atom optics 27 .

By introducing a more general mean-field which includes the GP one as a special case, we find another class of self-consistent excited states. Numerical examples show that the corresponding energies are substantially lower than those of the GP excited states. To avoid misunderstanding we mention that self-consistent excited states are physically distinct from collective (or particle-hole) excitations. The latter correspond to small oscillations around a given state and are described by mean-field linear-response theories based on the Bogoliubov approximation 28 ; 29 .

We consider a system of N identical bosons. The corresponding GP equation is the mean-field equation of this system with a $`\delta `$-function interaction potential $`W(\stackrel{}{r}_i\stackrel{}{r}_j)=\lambda _0\delta (\stackrel{}{r}_i\stackrel{}{r}_j)`$, where $`\stackrel{}{r}_i`$ is the position of the i-th boson and the nonlinear parameter $`\lambda _0`$ is related to the s-wave scattering length of the bosons 4 . The GP ansatz further assumes the wave function $`\mathrm{\Psi }_{GP}`$ to be a product of identical spatial one-particle functions (orbitals) $`\phi `$: $`\mathrm{\Psi }_{GP}(\stackrel{}{r}_1,\stackrel{}{r}_2,\mathrm{},\stackrel{}{r}_N)=\phi (\stackrel{}{r}_1)\phi (\stackrel{}{r}_2)\mathrm{}\phi (\stackrel{}{r}_N)`$. With this wavefunction the total energy $`E_{GP}=<\mathrm{\Psi }_{GP}|\widehat{H}|\mathrm{\Psi }_{GP}>`$ is obtained as the expectation value of our Hamiltonian $`\widehat{H}`$ and reads

$$E_{GP}=N\{\phi ^{}h\phi 𝑑\stackrel{}{r}+\frac{\lambda _0(N1)}{2}|\phi |^4𝑑\stackrel{}{r}\},$$
(1)
where $`h(\stackrel{}{r})=\widehat{T}+\widehat{V}(\stackrel{}{r})`$ is the unperturbed one-particle Hamiltonian consisting of the trap potential $`\widehat{V}(\stackrel{}{r})`$ and the kinetic energy operator $`\widehat{T}`$. By minimizing the energy (1) one obtains the GP equation:

$$\{h(\stackrel{}{r})+\lambda _0(N1)|\phi (\stackrel{}{r})|^2\}\phi (\stackrel{}{r})=\mu _{GP}\phi (\stackrel{}{r}).$$
(2)
Each self-consistent solution of this eigenvalue equation determines $`\phi `$ and the chemical potential $`\mu _{GP}`$. The ground state of the system, of course, corresponds to the solution with the lowest energy $`E_{GP}`$. The other solutions describe the self-consistent excited states of the system.

In the absence of the mutual interaction between the bosons, all bosons reside in a single spatial orbital in the ground state of the system. The GP ansatz for the wavefunction $`\mathrm{\Psi }_{GP}`$ is thus very appealing for the ground state and has indeed been very successful in explaining many observations. We have reason to assume, however, that not all relevant macroscopic excited states are describable by $`\mathrm{\Psi }_{GP}`$. The general wavefunctions $`\mathrm{\Psi }`$ of N non-interacting bosons is a product of orbitals which can all be different. Since the bosons are identical, this product must be symmetrized and several bosons may reside in the same spatial orbital. We may put $`n_1`$ bosons in orbital $`\varphi _1`$, $`n_2`$ in orbital $`\varphi _2`$ and so on. For transparency of presentation we restrict ourselves in the following to two orbitals $`\varphi _1`$ and $`\varphi _2`$ with particle occupations $`n_1`$ and $`n_2`$ and $`n_1+n_2=N`$. The extension to more orbitals is straightforward. The wavefunction now reads

$`\mathrm{\Psi }(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)=\widehat{𝒮}\varphi _1(\stackrel{}{r}_1)\mathrm{}\varphi _1(\stackrel{}{r}_{n_1})\varphi _2(\stackrel{}{r}_{n_1+1})\mathrm{}\varphi _2(\stackrel{}{r}_{n_1+n_2}),`$ (3)
where $`\widehat{𝒮}`$ is the symmetrizing operator.

By definition the most general mean-field energy is $`E=<\mathrm{\Psi }|\widehat{H}|\mathrm{\Psi }>`$ where $`\mathrm{\Psi }`$ is the above discussed wavefunction. It is evaluated to give 30

$`E=n_1h_{11}+n_2h_{22}+\lambda _0{\displaystyle \frac{n_1(n_11)}{2}}{\displaystyle |\varphi _1|^4𝑑\stackrel{}{r}}+`$
$`\lambda _0{\displaystyle \frac{n_2(n_21)}{2}}{\displaystyle |\varphi _2|^4𝑑\stackrel{}{r}}+2\lambda _0n_1n_2{\displaystyle |\varphi _1|^2|\varphi _2|^2𝑑\stackrel{}{r}},`$ (4)
where $`h_{ii}=\varphi _i^{}h\varphi _i𝑑\stackrel{}{r}`$ is the usual one-particle energy. We now minimize this energy with respect to the orbitals under the constraints that they are orthogonal and normalized, i.e.,

$`<\varphi _1|\varphi _2>=0,<\varphi _1|\varphi _1>=<\varphi _2|\varphi _2>=1.`$ (5)
This leads to the following set of coupled equations for the optimal orbitals:

$`\{h(\stackrel{}{r})+\lambda _0(n_11)|\varphi _1(\stackrel{}{r})|^2+2\lambda _0n_2|\varphi _2(\stackrel{}{r})|^2\}\varphi _1(\stackrel{}{r})=`$
$`=\mu _{11}\varphi _1(\stackrel{}{r})+\mu _{12}\varphi _2(\stackrel{}{r})`$
$`\{h(\stackrel{}{r})+\lambda _0(n_21)|\varphi _2(\stackrel{}{r})|^2+2\lambda _0n_1|\varphi _1(\stackrel{}{r})|^2\}\varphi _2(\stackrel{}{r})=`$
$`=\mu _{22}\varphi _2(\stackrel{}{r})+\mu _{21}\varphi _1(\stackrel{}{r}).`$ (6)
These equations should not be confused with those for condensates made of two types of bosons. The $`\mu _{ij}`$ are the Lagrange parameters due to the above mentioned constraints; $`\mu _{11}`$ and $`\mu _{22}`$ are related to the normalization constraints and the introduction of $`\mu _{12}`$ and $`\mu _{21}`$ enforces orthogonality of $`\varphi _1`$ and $`\varphi _2`$. We mention that $`n_1\mu _{12}=n_2\mu _{21}`$. In general, $`\mu _{12}0`$ unless $`\varphi _1`$ and $`\varphi _2`$ are of different spatial symmetry.

By setting $`n_2=0`$ we readily see that the general mean-field wavefunction $`\mathrm{\Psi }`$ in (3) contains the GP wavefunctions $`\mathrm{\Psi }_{GP}`$ as a special case. The energy $`E`$ in (4) then reduces to the expression for $`E_{GP}`$ in (1). Because of $`n_1\mu _{12}=n_2\mu _{21}`$, the parameter $`\mu _{12}`$ vanishes and the first equation in (6) is then nothing but the GP equation (2). Contrary to the GP equation (2), the general mean-field equations (6) are - owing to the presence of $`\mu _{12}`$ and $`\mu _{21}`$ \- not eigenvalue equations. This is a generic property of the latter equations. Unless either one of the particle occupations vanishes, or $`n_1=n_2=N/2`$, or $`\varphi _1`$ and $`\varphi _2`$ have different symmetry, or $`\mu _{11}=\mu _{22}`$, the off-diagonal Lagrange parameters $`\mu _{12}`$ and $`\mu _{21}`$ cannot be removed from (6) by linear transformations.

While we cannot attribute a physical observable to the off-diagonal Lagrange parameters $`\mu _{12}`$ and $`\mu _{21}`$ we can do so for the diagonal ones. As can be seen from (1), the energy needed to remove a boson from the condensate without changing the orbital $`\phi `$ is, within GP mean-field, given by $`E_{GP}(N)E_{GP}(N1)=\mu _{GP}`$. Analogously, we can compute within the present general mean-field the energy needed to remove a boson from orbital $`\varphi _1`$ and that from orbital $`\varphi _2`$. Recalling that $`\varphi _1`$ and $`\varphi _2`$ are orthogonal, we readily find from (4) and (6) the identities $`\mu _{11}=E(n_1,n_2)E(n_11,n_2)`$ and $`\mu _{22}=E(n_1,n_2)E(n_1,n_21)`$. Clearly, $`\mu _{11}`$ and $`\mu _{22}`$ can be viewed as chemical potentials of the $`\{\varphi _1\}`$\- and $`\{\varphi _2\}`$-boson manifolds, respectively.

Next, we can make use of the fact that the energy $`E`$ and the orbitals $`\varphi _1`$ and $`\varphi _2`$ depend on the particle occupation $`n_1`$ and treat this occupation as a variational parameter, i.e., search for its optimal value which makes the energy stationary. Interestingly, one can show for macroscopic occupancies $`n_1,n_21`$ that at any extremum of the energy as a function of $`n_1`$ the two quantities $`\mu _{11}`$ and $`\mu _{22}`$ coincide 30 :

$`\mu _{11}=\mu _{22}.`$ (7)
This is a relevant finding for understanding the concept of a condensate and the meaning of its self-consistent states. The system described by the wavefunction (3) consists of two (or more) subsystems each possessing its own chemical potential. This may contradict the picture one usually has of a condensate. It is thus important to note that these, in general different, chemical potentials become identical at the optimal occupations restoring thereby the picture of a condensate. With this in mind, we define as self-consistent condensate states of the present mean-field only those solutions of (6) at which the chemical potentials $`\mu _{11}`$ and $`\mu _{22}`$ are identical. These correspond to the extrema of the energy $`E`$ as a function of the boson occupation $`n_1`$.

Let us refer here to the problem of fragmentation often discussed in the literature, see refs. 31 ; 32  and references therein. A condensate is fragmented if its reduced one-body density matrix has two or more macroscopic eigenvalues (the total number of particles N is assumed large). Until now fragmentation has not been found for trapped condensates. Our present ansatz has the potential to describe fragmentation on the mean-field level. If the optimal boson occupations $`n_1`$ and $`n_2`$ are both macroscopic for large N, fragmentation in the respective self-consistent state indeed takes place. Consequently, we call such states which are beyond reach of the GP equation (2), self-consistent fragmented states.

We shall demonstrate in the following that the present mean-field supports self-consistent fragmented excited states. Moreover, these states can be at much lower energy than all the GP self-consistent excited states. As an example we choose a repulsive condensate $`(\lambda _0>0)`$ in an one-dimensional double-well potential. Since equations (6) have very recently been derived, the numerical procedures to evaluate them are not yet sufficiently developed to allow for computations in three dimensions. We have reason to believe, however, that similar results will be obtained in more dimensions.

In our example the potential

$$V(x)=c+\frac{\omega }{2}(x^2\sqrt{a^2+(\mathrm{\Delta }2xx_0)^2})$$
(8)
describes two wells separated by a barrier. The depths of these wells differ by a bias $`\mathrm{\Delta }`$. The parameter $`c`$ is chosen to set the bottom of the potential equal to zero. The coordinate $`x`$ is dimensionless and all energies and $`\lambda _0`$ are now in units of the frequency $`\omega `$. The equations for the optimal orbitals have been evaluated numerically by using the DVR method 33 . Starting from an initial guess for $`\varphi _1`$ and $`\varphi _2`$, for instance, the solutions of the equations with slightly different boson occupations, the equations are iterated until self-consistency is achieved. To accelerate the convergence of the calculations we have employed an energy-shift technique as often done in Hartree-Fock electronic structure computations 34 . The GP equation has been evaluated by solving the corresponding eigenvalue equation self-consistently using the DVR method.

Having obtained the optimal orbitals for the present and GP mean-fields, we have computed the energies $`E`$ and $`E_{GP}`$. Computations have been carried out systematically for the ground and the first few excited states using various values of the coupling constant $`(\lambda =\lambda _0N)`$ (note that $`N1`$) and of the parameters appearing in the potential $`V(x)`$ in (8). For each set of these values the optimal orbitals and the mean-field energy $`E`$ have been determined as a function of the boson occupation $`n_1`$. A typical example is shown in fig.1. In this figure the mean-field energy per particle $`E/N`$ is depicted as a function of the fractional occupation $`n_1/N`$ for $`\lambda =2.5`$. For comparison, the energies per particle $`E_{GP}/N`$ of the GP ground and first excited state are indicated as well. Two disjoint $`E/N`$ curves are seen in the figure and we concentrate first on the one at lower energy. Starting the calculations at $`n_1=N`$, the energies $`E`$ and $`E_{GP}`$, of course, coincide. We then continuously lowered the value of $`n_1/N`$ and obtained a smooth $`E/N`$ curve. This curve exhibits two maxima and a minimum which is, of course, located between them. Following (7) there is only one value for the chemical potential at each of these extrema. Accordingly, the states corresponding to these extrema are self-consistent fragmented excited states of the condensate.

The two maxima of $`E(n_1)`$ correspond to metastable states. Any variation of the boson occupation $`n_1`$ lowers the energy favoring a ”decay” into either the GP ground state which is also the ground state of the present mean-field or to the state corresponding to the minimum of $`E(n_1)`$. The latter state is stable with respect to variations of $`n_1`$, or briefly a stable state. A substantial energy barrier has to be overcome in order to lower the energy and to destabilize this state. Using equations (4-6) we can show that

$$\frac{dE}{dn_1}=\mu _{11}\mu _{22}$$
(9)
holds. This physically appealing relation adds insight into the understanding of the curves $`E/N`$ depicted in fig.1. For large N the quantity $`\mu _{11}\mu _{22}`$ is the energy gain obtained by removing a boson from the $`\varphi _1`$-manifold and adding it to the $`\varphi _2`$-manifold. $`\mu _{11}\mu _{22}`$ is thus the ”driving force” for the flow of bosons between the two boson subsystems. At the extrema of $`E`$ as a function of the boson occupation one recovers (7).

The orbitals $`\varphi _1(x)`$ and $`\varphi _2(x)`$ for the stable fragmented excited state discussed above are depicted in fig.2A. For comparison also the orbitals $`\phi `$ for the ground and for the first excited GP state are shown. The ground state GP orbital is as usual nodeless and consists of two humps related to the two wells of the potential. The first excited state GP orbital exhibits a node close to the top of the barrier. Each of the orbitals $`\varphi _1`$ and $`\varphi _2`$ of the fragmented stable excited state also possesses a node, but these nodes are located at distinct different sites. $`\varphi _1`$ is rather localized at the deeper well and its node is close to the minimum of the other well. The reverse situation holds for $`\varphi _2`$. That these nodes are well separated can be understood by comparing the density per particles $`\rho =(n_1|\varphi _1|^2+n_2|\varphi _2|^2)/N`$ with the respective quantities $`\rho _{GP}=|\phi |^2`$ for the GP states. These densities are shown in fig.2B together with the trap potential $`V(x)`$. Since the bosons in our example repeal each other, the density would like to spread over space to reduce this repulsion. Indeed, in the ground state the density is substantial in between the wells. Due to the node of the excited GP state, the density vanishes close to the top of the barrier, and as a consequence the density is enhanced inside the wells. This inevitably leads to a substantial increase of the energy $`E_{GP}`$ for the excited state. Fragmentation assigns different nodes to different orbitals and hence the density can penetrate the region of the barrier as seen in fig.2B. Consequently, the energy of the fragmented excited state is much lower than that of the GP excited state as seen in fig.1.

It should be stressed that the above findings depend only weakly on the shape of the double-well potential. We have performed numerous calculations varying the coupling constant $`\lambda `$, the bias $`\mathrm{\Delta }`$ and other parameters of the trap potential. In particular, if we choose the potential to be symmetric (bias $`\mathrm{\Delta }=0`$), the same conclusions as above can be drawn. Furthermore, enlarging the coupling constant $`\lambda =\lambda _0N`$ results in even more prominent minimum and maxima of the energy curves as a function of $`n_1/N`$.

Let us briefly discuss the short $`E/N`$ curve in fig.1. We have obtained this curve by solving (6) starting from the GP excited state orbital, i.e., from $`n_2=N`$, and increasing $`n_1`$ continuously. As seen in the insert which shows this curve on an enlarged scale, the curve exhibits a maximum thus describing another fragmented metastable state. It is also seen that the curve has a boundary minimum at $`n_1=0`$, indicating that the corresponding GP excited state is stable. The energy barrier involved is, however, very low. We remark that at large values of $`\lambda (8.0)`$, this barrier disappears and we find a boundary maximum, i.e., the GP state becomes metastable. The $`E(n_1)`$ curve then acquires a minimum at $`n_10`$ giving rise to a new stable fragmented state.

In conclusion, self-consistent fragmented excited states of condensates exist. These states can be either metastable or stable with respect to the flow of bosons from one fragment (manifold of bosons residing in one orbital) to another one. In the class of examples studied here energies of fragmented states are much lower than those of the excited GP states. Due to the transparent physics behind populating several orbitals in excited states, we expect such states to exist also in two and three dimensions and in other trap potentials.

The authors acknowledge useful discussions with Ofir Alon.


format: text

text:
# Enhancement of the escape time in metastable states with colored noise

## 1 Introduction

Nonlinear relaxation decay of physical systems from an initial unstable or metastable state involves fundamental aspects of non-equilibrium statistical mechanics. The investigation of nonlinear dynamics and instabilities in systems away from equilibrium in fact has led to some counterintuitive and resonance-like phenomena. Among these we cite the stochastic resonance , the resonant activation , the noise enhanced stability (NES)  and noise induced phase transitions . Recent theoeretical investigations have shown that the average escape time from metastable states in fluctuating potentials has a nonmonotonic behaviour as a function of the noise intensity . This resonance-like behaviour, which contradicts the monotonic behaviour predicted by the Kramers theory , is the NES phenomenon: the stability of metastable or unstable states can be enhanced by the noise and the average life time of the metastable state is greater than the deterministic decay time.

The inclusion of realistic noise sources with a finite correlation time in modelling dynamical systems can impact both the stationary and the dynamic features of nonlinear systems. For metastable thermal equilibrium systems it has been demonstrated that colored thermal noise can substantially modify the diffusive barrier transmission . For bistable systems the colored noise driven escape rate for small and large correlation time has been derived in ref. . A rich and enormous literature on escape processes driven by colored noise was produced in the 80’s. As an example of this literature concerning colored noise we cite the escape from metastable states , the MFPT from a marginal state  and the decay of an unstable state . In this work we present a study of the average decay time of an overdamped Brownian particle subject to a cubic potential with a metastable state. We focus on the role of different unstable initial conditions and of colored noise in the average escape time. The effect of *color* on the transient dynamics of the escape process is strictly related to the characteristic time scale of the system, i. e. the relaxation time inside the metastable state $`\tau _r`$. For correlation time values less than the relaxation time ($`\tau _c<\tau _r`$), the system ”sees” the noise source as white noise and the dynamical regime of the Brownian particle is like the white noise dynamics with a shift of all curves calculated at different initial positions towards higher values of noise intensity. For strong color, i. e. for values of correlation time larger than the relaxation time ($`\tau _c>\tau _r`$) nontrivial results are obtained. Particularly we obtain: (i) a big shift of the behaviours of the average escape times towards higher noise intensities; (ii) an enhancement of the value of the average escape time maximum and a broadening of the NES region in the plane ($`\tau ,D`$), which becomes very large for high values of $`\tau _c`$; (iii) the absence of the peculiar initial position $`x_c`$, found in a previous study , which separates the set of the initial unstable states producing divergency from those which give only a nonmonotonic behavior of the average escape time. This is the relevant result of this paper. For all the initial unstable states we find an enhancement of the average escape time with respect to the deterministic decay time and in the strong correlated noise regime ($`\tau _c>>\tau _r`$) a nonmonotonic behavior of the the average escape time.

## 2 The model

The starting point of our study is the Langevin equation

$$\dot{x}=\frac{U(x)}{x}+\eta (t)$$
(1)
where $`\eta (t)`$ is the Ornstein-Uhlenbeck process

$$d\eta =k\eta dt+k\sqrt{D}dW(t)$$
(2)
where $`W(t)`$ is the Wigner process giving a Gaussian noise with the usual statistical properties: $`\xi (t)=0`$ and $`\xi (t)\xi (t+\tau )=\delta (\tau )`$. The integration of Eq.(2) yields

$$\eta (t)=\eta (0)e^{kt}+k\sqrt{D}_0^te^{k(tt^{})}𝑑W(t^{})$$
(3)
which for $`\eta (0)=0`$ and for $`t\mathrm{}`$ gives

$$\eta (t)=0,$$
(4)
$$\eta (t)\eta (t+\tau )=\frac{kD}{2}e^{k\tau },$$
(5)
where $`1/k=\tau _c`$ is the correlation time of the process. The system of stochastic differential equations (1) and (2), which represent a two-dimensional Markovian process, is equivalent to a non-Markovian Langevin equation driven with additive Gaussian correlated noise with $`\eta (t)`$ obeying the properties (4) and (5).

It is possible to show that the Eq.(2) for the OU process gives the correct limit $`lim_{\tau _c0}\eta (t)=\sqrt{D}\xi (t)`$, i. e. the white noise term. In fact

$`\underset{\tau _c0}{lim}\eta (t)`$ $`=`$ $`2\sqrt{D}{\displaystyle _0^t}\underset{\tau _c0}{lim}{\displaystyle \frac{e^{(tt^{})/\tau _c}}{2\tau _c}}dW(t^{})`$ (6)
$`=`$ $`2\sqrt{D}{\displaystyle _0^t}\delta (tt^{})\xi (t^{})d(t^{})=\sqrt{D}\xi (t),`$
where we use the property of $`\delta `$ function

$$_{\mathrm{}}^+\mathrm{}𝑑x\delta (xx_o)f(x)=f(x_o),$$
(7)
and the factor $`2`$ disappears because of the integration range. The stationary correlation function of the Ornstein-Uhlenbeck process (5) gives in the limit $`\tau _c0`$ the correlation function of the white noise

$$\underset{\tau _c0}{lim}\eta (t)\eta (t+\tau )=\underset{\tau _c0}{lim}\frac{D}{2\tau _c}e^{\tau /\tau _c}=D\delta (\tau ).$$
(8)
The potential $`U(x)`$ used in Eq.(1) is a cubic one

$$U(x)=ax^2bx^3,$$
(9)
with $`a=0.3`$, $`b=0.2`$. The potential profile has a local stable state at $`x_0=0`$ and an unstable state at $`x_0=1`$ (see Fig.1).

The relaxation time for the metastable state at $`x_0=0`$ is

$$\tau _r=\left[\frac{d^2U(x)}{dx^2}\right]_{x=0}=2a,$$
(10)
which is the characteristic time scale of our system.

## 3 Results and Comments

The calculations of the average escape time as a function of the colored noise intensity have been performed by averaging over $`5000`$ realizations the stochastic differential equation (1). The absorbing boundary for the escape process is put on $`x_F=20`$. For all the initial unstable states (see Fig.1) and all the correlation times evaluated we find an enhancement of the average escape time with respect to the deterministic time as a function of the noise intensity. We can observe the characteristic behavior of the NES effect for unstable states with white noise in Fig.2a  where we also find the comparison with the calculation performed with colored noise with $`\tau _c=0.01`$ (Fig.2b). The behaviours of the average escape time as a function of colored noise intensity with the others values of $`\tau _c`$ are shown in Fig.3.

We clearly observe two dynamical regimes: (a) weak colored noise ($`0<\tau _c<\tau _r`$) and (b) strong colored noise ($`\tau _c>\tau _r`$). Starting our observation from Fig.3 we see that the qualitative behavior of MFPT with white noise is recovered with $`\tau _c=0.1`$, until $`\tau _c\tau _r=0.6`$. That is an enhancement of the average escape time with respect to the deterministic decay time with a divergent behaviour for $`x_{max}<x_0<x_c`$ and a non monotonic behavior for $`x_0x_c`$.

Increasing the value of the correlation time Fig.3 ($`\tau _c\tau _r`$) we observe a shift of the maximum towards higher values of the noise intensity and an enhancement of the value of the average escape maximum. Moreover we observe a broadening of the NES region, which becomes very large (up to twelve order of magnitude) for high values of the correlation time $`\tau _c`$ (see Fig.3 with $`\tau _c=100`$). The unexpected result is that increasing the correlation time the divergent behaviour tends to disappear, as shown in Fig.3 for $`\tau _c=1.0`$ and initial position $`x_o=1.73`$, and disappear for all the initial positions investigated for $`\tau _c=10`$ and $`\tau _c=100`$.

For high values of the noise intensity all the graphs show a monotonic decrease behavior as a function of noise intensity collapsing in a unique curve, according to the behavior predicted by the Kramers theory . Moreover the slope of this limit curve becomes flatter increasing the correlation time. This means that the NES effect involves more and more orders of magnitude of the noise intensity. The effect of the colored noise is therefore to delay the escape process or in other words to enhance more and more the stability of the metastable state.

## 4 Conclusions

In this work we analyzed the effect of the colored noise, generated by an OU process, on the enhancement of the average escape time in a cubic potential with a metatstable state. We analyze different initial unstable states with $`x_0>x_c`$, where $`x_c`$ is the crossover between the divergent behaviour and the nonmonotonic behaviour of the average escape time . We obtain NES effect for all the initial positions investigated and an enhancement of the NES effect for increasing values of correlation times. The results obtained for a particle moving in a cubic potential are quite general, because we always obtain NES effect when a particle is initially located just to the right of a local potential maximum, with a local minimum or metastable state in its left side and the global escape region in its right side.

From a comparison with analogue calculation with white noise, we observe the following similarity and differences for increasing values of the correlation times: 1) the peculiar point $`x_c`$ found for white noise separating divergent behaviours from non-monotonic behaviour, is present only in the weak colored regime $`\tau _c<\tau _r`$; 2) all the divergencies present in the case of white noise for all the unstable initial position in the range $`x_{max}<x_0<x_c`$ tend to transform in pure NES effect, i. e. without divergencies, for increasing values of $`\tau _c`$; 3) the region in the plane $`(\tau ,D)`$ of the NES effect increases up to $`12`$ order of magnitude for increasing values of $`\tau _c`$ (strong colored noise regime).

Nonmonotonic behavior of the mean escape time as a function of noise intensity therefore is a noise-induced effect for nonlinear nonequilibrium systems with metastable states which is enhanced when we consider realistic noise sources with finite correlation time.

In experiments real noise sources are correlated with a finite correlation time. As a consequence the NES effect can be observed at higher noise intensities with respect to the idealized white noise case. Therefore the theoretical results obtained using white noise source are only a good approximation for low values of the correlation time. The enhancement and the shift of the NES region towards higher values of the noise intensity allow to reveal experimentally the NES effect using a suitable correlation time $`\tau _c`$.

## 5 Acknowledgments

This work was supported by INTAS Grant 01-0450, by INFM and MIUR.


text:
# Non-Gaussianity of Resistance Fluctuations Near Electrical Breakdown

## I Introduction and Model

Electrical breakdown is an irreversible change in the resistance of a two terminal device and it usually occurs in the presence of high applied electric stress, like current or voltage. Breakdown of condunctor-insulator composites , soft dielectric breakdown in ultrathin oxide films  and degradation of metallic films due to electromigration  are typical experimental examples. Nonlinear current voltage characteristics, giant enhancement of excess noise over the quadratic dependence on the applied field are the mostly studied breakdown precursors. . Yet another important and less understood breakdown precursor is represented by the emergence of non-Gaussian fluctuations , which origin has attracted an increasing interest in the recent literature . For these reasons, here we study the distribution of the resistance fluctuations of a thin film near electrical breakdown. To this purpose we make use of the Stationary Network Under Biased Percolation (SNUBP) model . This model allows the study of the electrical conduction of disordered materials over the full range of applied bias values, from the linear regime up to the breakdown. Moreover, the SNUBP provides a good modeling of many features associated with the electrical instability of composites materials  and with the electromigration damage of metal lines , two important classes of breakdown phenomena.

In the SNUBP model , the resistance and the resistance fluctuations of a disordered film are determined by the competition of two biased processes taking place in a two-dimensional square-lattice resistor network. More precisely, the network consists of $`2N^2`$ resistors in two possible states: (i) regular, corresponding to resistors with resistance $`r_n=r_0[1+\alpha (T_nT_0)]`$ and (ii) broken, corresponding to resistors with $`r_{OP}=10^9r_n`$. Here, $`N`$ determines the linear sizes of the network, $`\alpha `$ is the temperature resistance coefficient, $`T_0`$ the bath temperature, and $`T_n`$ the local temperature of the n-resistor, resulting from Joule heating and thermal exchange with neighbour resistors :

$$T_n=T_0+A\left[r_ni_n^2+\frac{3}{4N_{neig}}\underset{l=1}{\overset{N_{neig}}{}}\left(r_li_l^2r_ni_n^2\right)\right]$$
(1)
where $`N_{neig}`$ is the number of first neighbors around the nth resistor and $`i_n`$ the current flowing in it. The parameter $`A`$ represents the thermal resistance of each resistor and it determines the importance of Joule heating effects. The network is contacted at the left and right hand sides to perfectly conducting bars through which a constant stress current $`I`$ is applied. The two biased processes consist of stochastic transitions between the two possible states and they occur with thermally activated probability, $`W_D`$ and $`W_R`$, which depend on the local temperature and are characterized by two activation energies, $`E_D`$ and $`E_R`$, with $`E_D>E_R`$ for networks of finite sizes . The network evolution is obtained by a Monte Carlo simulation which updates the network resistance after breaking and recovery processes according to an iterative procedure as detailed in Ref. . The sequence of the successive configurations provides the final $`R(t)`$ signal with the time scale calibrated on the iteration steps. Then, depending on the $`E_D`$ and $`E_R`$ values and on the stress conditions, the network reaches a steady state or undergoes an irreversible electrical failure. This last possibility is associated with the existence of at least one continuous path of broken resistors between the upper and lower sides of the network, i.e. with the achievement of the percolation threshold, $`p_C`$, for the fraction of broken resistors. It has been found  that for a given value of $`T_0`$, a threshold current value, $`I_B`$, exists above which the electrical breakdown occurs. On the other hand, for $`II_B`$ the network is stable, i.e. the fraction of broken resistors, $`p`$, and the network resistance, $`R`$, fluctuate around their average values $`<p>`$ and $`<R>`$, respectively. In the following section we will analyze the results of simulations carried out by considering square networks of size $`75\times 75`$ and by taking the following values for the other parameters (chosen as physically plausible): $`r_0=1`$ ($`\mathrm{\Omega }`$), $`\alpha =10^3`$ (K<sup>-1</sup>), $`A=5\times 10^5`$ (K/W), $`E_D=0.170`$ (eV), $`E_R`$ in the range $`0.026÷0.155`$ (eV) and $`T_0=300`$ (K). The values of the external current range from $`0.001I3.0`$ (A).

## II Results

Figure 1 reports a set of resistance evolutions for increasing current values. In this case the activation energy for the recovery process is $`E_R=0.043`$ (eV). This value of $`E_R`$ implies a breakdown current $`I_B=2.1`$ (A). The figure displays two important features of the electrical response of conducting films to high fields: i) the average resistance exhibits an Omic behavior at low currents and then increases significantly at high currents (non-Ohmic regime); ii) the amplitude of resistance fluctuations increases strongly for currents near to the breakdown value. By analyzing the first and the second moment of the steady resistance we have found that the non ohmic regime can be divided in two regions . Actually, a moderate bias region, where both the average resistance $`<R>`$ and the relative variance of resistance fluctuations $`\mathrm{\Sigma }<(\mathrm{\Delta }R)^2>/<R>^2`$ increase quadratically with the current, is followed by a pre-breakdown region where superquadratic behaviors emerge .

Remarkably, another feature of the resistance fluctuation distribution is found to emerge: the non-Gaussian behavior which accounts for the importance of moments higher than the second. Here, we have thus investigated the Gaussianity of the $`R(t)`$ steady-state signals shown in Fig. 1. We have considered time series containing about $`3\times 10^5`$ resistance values. The results of the analysis are reported in Fig. 2 which shows the probability density function (PDF), $`\varphi `$, of the distribution of resistance fluctuations for two values of the bias current: $`I=0.01`$ (A) (triangles) and $`I=I_B=2.1`$ (A) (full circles). Precisely, in this figure we have reported on a lin-log plot the product $`\sigma \mathrm{\Phi }`$ as a function of $`(<R>R)/\sigma `$ (where $`\sigma `$ is the root mean square deviation). In fact, by making the distribution independent of its first and second moments, this normalized representation is particularly convenient to explore the functional form .

For comparison, in Fig. 2 we also show the Gaussian distribution (long dashed line), which in this normalized representation has zero mean and unit variance. We notice that in 1998 Bramwell, Holdsworth and Pinton  realized that the PDF of a global quantity of a system at a critical point takes a universal behavior, irrespectively of the particular quantity and of the nature and sizes of the system. Successively, the functional form describing this universal behavior was identified . Therefore, in Fig. 2 we report also the Bramwell, Holdsworth and Pinton (BHP) distribution (thick continuos line). Before giving its expression, we have to introduce the following definitions . We call $`m`$ a fluctuating quantity (for example the magnetization of a ferromagnet like in Ref. ), $`<m>`$ and $`\sigma _m`$ its mean value and root mean square deviation respectively, $`P(m)`$ its PDF, $`y(m<m>)/\sigma _m`$ the normalized variable, $`\mathrm{\Pi }(y)\sigma _mP(y)`$ the normalized PDF and finally $`xb(ys)`$. Then, the BHP distribution takes the following expression :

$$\mathrm{\Pi }(y)=K[e^{xe^x}]^a$$
(2)
where $`a=\pi /2`$, $`b=0.936\pm 0.002`$, $`s=0.374\pm 0.001`$ and $`K=2.15\pm 0.01`$. Equation (2) can be considered a generalization of the Gumbel distribution, which is often associated with the occurrence of rare events. Thus, looking at Fig. 2 we can outline the following main points: i) the resistance fluctuation distribution exhibits a non-Gaussian behavior; ii) the non-Gaussianity increases for currents near to the breakdown value; iii) even at low biases (i.e within the Ohmic regime) the distribution keeps a weak non-Gaussianity; iv) the distribution at $`I_B`$ departs significantly from the universal behavior. Point iii) is understood as a finite size effect, as confirmed by the analysis of data for networks of different sizes, not reported here. Thus, the points iii) and iv), by indicating a size dependent and non universal behavior, both witness a first order transition, in agreement with the results in Refs.  and with electrical breakdown experiments in composites . On the other hand, by recalling that real systems have finite sizes, the feature ii) can provide an important precursor of failure for systems of finite and given sizes. In order to identify what kind of distribution describes the data in Fig. 2, we have considered several distributions among those available in the literature. Accordingly, we have found that the best description to the present data is given by the expression in Eq. (2) once the parameters $`a`$, $`b`$, $`s`$ and $`K`$ are considered as fitting parameters. The dashed and the solid curves in Fig. 2 represent the best-fit with Eq. (2) to the data for $`I=2.1`$ (A) (full-circles) and for $`I=0.01`$ (A) (triangles), respectively.

In the following, we further consider the effect of the activation energy of the defect recovery process on the PDF. In fact, the recovery energy, by controlling the level of disorder inside the system, plays an important role on the system behavior . To this purpose, Fig. 3 reports the distributions of the resistance fluctuations obtained for different values of the recovery energy: $`E_R=0.043`$ (eV) (full circles), $`E_R=0.130`$ (eV) (stars), $`E_R=0.155`$ (eV) squares. In all the cases, the current corresponds to the breakdown, and the values are: $`I=2.1`$ (A), $`I=0.38`$ (A) and $`I=0.095`$ (A), respectively. Thus, full circles are the same data of Fig. 2. We notice that $`E_R=0.155`$ (eV) is very close to the $`E_D`$ value and practically represents one the largest value of $`E_R`$ yet providing a steady state of the network . Therefore, as discussed in Ref. , for this value of $`E_R`$ the network is very near to its critical point and it is characterized by an average fraction of defects ($`<p>=0.35`$) close to the value of the random percolation threshold (0.5 for bond percolation on a square and infinite lattice). Accordingly, from Fig. 3 we can see that the resistance fluctuation distribution corresponding to this value of $`E_R`$ is well described by the universal BHP distribution. The solid and the dashed curves in this figure represent the best-fit with Eq. (2) to the data obtained for $`E_R=0.130`$ (eV) and $`E_R=0.043`$ (eV), respectively. Also in these cases the fit with Eq. (2) is found to be satisfactory.

In conclusions, we have studied the distribution of resistance fluctuations of a thin film under different bias conditions and for several values of the recovery energy. The distribution exhibits a non-Gaussian behavior which is well described by a generalized form of the Gumbel distribution. Furthermore, a value of the recovery energy close to that of defect generation, a condition which corresponds to an average defect fraction approaching the random percolation value, is identified as the condition under which the distribution achieves the universal behavior described by the BHP distribution .

## ACKNOWLEDGEMENTS

This research has been performed within the SPOT-NOSED project IST-2001-38899 of the EC.


text:
# Comparative study of defect energetics in HfO2 and SiO2

## Abstract

We perform ab initio calculations, based on density functional theory, of substitutional and vacancy defects in the monoclinic hafnium oxide (m-HfO2) and $`\alpha `$-quartz (SiO2). The neutral oxygen vacancies and substitutional Si and Hf defects in HfO2 and SiO2, respectively, are investigated. Our calculations show that, for a large range of Hf chemical potential, Si substitutional defects are most likely to form in HfO2, leading to the formation of a silicate layer at the HfO2/Si interface. We also find that it is energetically more favorable to form oxygen vacancies in SiO2 than in HfO2, which implies that oxygen deficient HfO2 grown on top of SiO2 will consume oxygen from the SiO2.

The continuous device miniaturization in the microelectronic industry will eventually lead, within the present technology, to the end of the use of amorphous SiO2 (a-SiO2) as gate dielectric in metal-oxide-semiconductor field-effect transistors (MOSFETs). The existence of a thickness limit for the a-SiO2 around 10-12 Å, has clearly been established experimentally Muller . One way to circumvent this problem, still keeping Si as the basic device material, is to employ high-permitivity materials as alternative gate dielectrics in place of the conventional a-SiO2. Among them, hafnium oxide is emerging as the material with greatest potential to substitute SiO2, mainly due to its high dieletric constant and thermodynamic stability, when it forms interface with Si.

Even though hafnium oxide is thermodynamically stable against an overall decomposition as Hf and SiO2 when grown on Si, interfacial reactions can occur. Thereby the formation of a thin interfacial layer (oxides, silicates and silicides) between the HfO2 and the Si surface, has been recently observed Roh ; Jeon . This interfacial layer occurs during almost any film growth processes or post-annealing, which is an intrinsic part of any growth cycle. Therefore, the thermodynamic stability of the hafnium oxide in contact with silicon is identified as a critical issue for the application of alternative gate dielectric in silicon-based devices Misa ; Bush . Moreover, the study of possible defects related to the migration of atoms across the interface is of fundamental importance. In particular, a significant source of defects in this system is the interface itself, which has been shown Gutowski  to consist of Hf silicates with dielectric constant lower than that of HfO2 Callegari ; Kato .

In the present work we address the formation of neutral defects through first-principles calculations, based on the density functional theory (DFT). We analyze the formation of Si substitutional defects in HfO2, as well as Hf substitutional defects in SiO2, for different growth conditions. Finally, the energetics of an oxygen vacancy in SiO2 is compared to a similar vacancy in HfO2, in order to understand the growth of hafnium oxide under oxygen-poor conditions.

Many experimental works Cho ; Park ; Moonju ; Bastos  have addressed the chemical reactions that could occur in the HfO2/Si interface during the HfO2 growth cycle. This is due in order to prevent and/or control the interfacial layer formation. Almost all works have reported the formation of an interfacial Hf silicate in oxygen-rich atmospheres. Furthermore, Wang and co-workers Lim  have shown that in opposite conditions, i.e., oxygen-deficient atmospheres, the Hf silicate interfacial formation does not happen during the HfO2 growth cycle. However, there are still many open questions, such as which atomic species are migrating when the interfacial silicate is formed. For example, it is important to know if Hf will be incorporated in a formed SiO2 layer, or if Si from either the bulk Si or this SiO2 layer will be incorporated in the HfO2.

The DFT calculations were performed using ultrasoft Vanderbilt pseudopotentials Vanderbilt , and the generalized gradient approximation(GGA) for the exchange-correlation potential as implemented in VASP code Perdew ; vasp1 ; vasp2 ; vasp3 . In order to study the defects in the different systems we have considered the monoclinic HfO2 and the $`\alpha `$-quartz (SiO2) crystalline phases, using a 96 atoms and 72 atoms supercells, respectively. For these cells, we have used a plane wave cutoff energy of 400 eV and a 2$`\times `$2$`\times `$2 Monkhorst-Pack k-mesh. These crystal structures have been previously used to describe these systems nosso1 ; nosso2 ; foster . In all calculations the atoms were allowed to relax until all components of the atomic forces were smaller than 0.025 eV/Å.

The silicon substitutional defect (Si$`_{\text{Hf}}`$) was created in the m-HfO2 supercell by substitution of one hafnium atom by one silicon atom in the equilibrium perfect crystal. On the other hand, the hafnium substitutional defect (Hf$`_{\text{Si}}`$) was created in the $`\alpha `$-quartz by substitution of one silicon atom by one hafnium atom in the equilibrium perfect crystal. Our results show that the presence of either Si$`_{\text{Hf}}`$ or Hf$`_{\text{Si}}`$ in HfO2 or SiO2, respectively, do not introduce any additional active levels in the band gap. For the perfect HfO2 at the equilibrium, we obtain an indirect band gap of 3.9 eV along $`\mathrm{\Gamma }`$-B. For this same supercell, the density of states (DOS) shows that the O(2s) and O(2p) bands are centered at around -14 and 0 eV, with bandwidths of approximately 4 eV and 8 eV, respectively, below valence-band edge, whereas the Hf(5d) band is centered at around 8 eV forms in the conduction band, with bandwidth of 6 eV. These results are in good agreement with previous DFT calculations of this material kang . Moreover, for HfO2 with a Si substitutional defect, we observe a resonant state around 6.9 eV below the top of the valence band, which is related to the O(2p)-like and Si bonding state Binggeli .

The formation energy for a Si$`_{\text{Hf}}`$ is calculated as

$`E_f^{Si_{\text{Hf}}}=[E_t(Si_{\text{Hf}})+\mu _{Hf}][E_t(HfO_2)+\mu _{Si}],`$ (1)
whereas for Hf$`_{\text{Si}}`$ the similar expression is

$`E_f^{Hf_{\text{Si}}}=[E_t(Hf_{\text{Si}})+\mu _{Si}][E_t(SiO_2)+\mu _{Hf}].`$ (2)
In the above expressions, Et($`D_S`$) are the total energies of the fully relaxed supercells (either m-HfO2 or SiO2) with the substitutional defect $`D_S`$, and Et($`XO_2`$) are the total energies of the similar supercells for the perfect crystals of XO2 (X = Hf or Si). The values of the chemical potentials, $`\mu _{Hf}`$ and $`\mu _{Si}`$, depend on the growth conditions. We have considered two limits for $`\mu _{Hf}`$: (i) The bulk metal as a reference, which would correspond to a Hf-rich growth condition and the formation of Hf clusters in the bulk or at the surface of the oxide; (ii) Under oxygen rich conditions, or considering that there is plenty of oxygen atoms in the HfO2 nieminen , such that the removed hafnium remains always in equilibrium with the oxygen gaseous, the $`\mu _{Hf}`$(gas) can be obtained as $`\mu _{Hf}`$(gas) = $`\mu _{HfO_2}\mu _{O_2}`$. The HfO2 chemical potential ($`\mu _{HfO_2}`$) was obtained as the energy per unit formula for the monoclinic bulk hafnium oxide, and $`\mu _{O_2}`$ is the energy of an isolated oxygen molecule, which was obtained through a DFT total energy calculation for an O2 inside a cubic supercell of 15 Å side. For the silicon chemical potential, two similar limits were considered: (i) $`\mu _{Si}`$ as the crystalline bulk Si chemical potential, and (ii) $`\mu _{Si}`$ as the chemical potential of Si in SiO2 under oxygen rich conditions, i.e., $`\mu _{Si}`$ = $`\mu _{SiO_2}\mu _{O_2}`$, with $`\mu _{O_2}`$ as above and $`\mu _{SiO_2}`$ obtained as the energy per unit formula of $`\alpha `$-quartz.

The results for the formation energies are presented in Fig. 1. As can be seen from Fig. 1(a), under oxygen rich conditions, the formation of a Si substitutional defect in HfO2 is very likely to occur, specially if the Si source is the bulk crystal. Under this O-rich situation, the formation of this defect would only become unfavorable if a thick enough layer of SiO2 would exists between the Si substrate and the HfO2. In this Case, either the source of Si atoms would be the silica, or the Si would somehow not diffuse through this thick layer. If the SiO2 is thin enough, such that the Si diffusion through it is significant, then the formation of a Hf silicate (at least close to the interface) seems unavoidable lee . On the other hand, under Hf rich conditions, it is less likely that Si substitutional defects in HfO2 will be formed.

The formation of Hf substitutional defects in SiO2 has a behavior that is opposite to Si$`_{\text{Hf}}`$, as can be seen in Fig. 1(b), i.e., it is not likely to form under oxygen rich conditions, and becomes more probable under Hf-rich conditions. In a situation where the Si chemical potential is given by its bulk value (either thin or no SiO2 layer), the formation of a Hf$`_{\text{Si}}`$ in SiO2 is unlikely to happen, for a large range of Hf chemical potential.

We also considered the formation of neutral oxygen vacancies ($`V_O`$), both in SiO2 as well as in HfO2 foster , since they can be created in films and bulk samples due to the growth cycle. The neutral oxygen vacancy in the m-HfO2 and $`\alpha `$-quartz were generated by simple removal of an oxygen atom, followed by full relaxation of all remaining atoms. The formation energies for a $`V_O`$ in XO2 (X=Hf or Si), $`E_f(V_O)`$, were calculated as,

$`E_f^{XO_2}(V_O)=[E_t^{XO_2}(V_O)+\mu _O][E_t(XO_2)],`$ (3)
where $`E_t^{XO_2}(V_O)`$ and $`E_t(XO_2)`$, are the total energies of supercells of $`XO_2`$ (X = Hf or Si) with, and without, an oxygen vacancy, respectively. The oxygen chemical potential, $`\mu _O`$, was considered either as the total energy of an isolated oxygen atom, or as one half of the energy of an isolated oxygen molecule (in both cases, the O2 triplet ground state was used). The monoclinic phase has non-equivalent oxygen atoms, i.e, which some sites are threefold coordinated, whereas others are fourfold coordinated by hafnium atoms. In this way, we have determined the formation energy for both vacancies types using Eq. (3) and we obtain a formation energy difference around 0.02 eV. In the former case, we obtained of $`E_f^{HfO_2}(V_O)`$ = 9.32 eV and $`E_f^{SiO_2}(V_O)`$ = 8.10 eV, whereas for the latter choice of $`\mu _O`$, we obtained $`E_f^{HfO_2}(V_O)`$ = 6.38 eV and $`E_f^{SiO_2}(V_O)`$ = 5.16 eV. This indicates that, although neutral oxygen vacancies are energetically unfavorable in both materials, they are more stable in silicon oxide than in hafnium oxide, by approximately 1.23 eV. This implies that, if an oxygen deficient HfO2 is grown on top of a SiO2 layer, oxygen atoms will migrate from SiO2. In this way, oxygen vacancies are created in the silica layer toward the HfO2, healing at the sometimes oxygen vacancies in the hafnia. Indications that this process does indeed occur, has been recently reported Lim .

In summary, our results show that, unless the hafnium chemical potential is always very close to its bulk value,i.e., oxygen-poor growth conditions, the formation of Si substitutional defects in the HfO2 is almost unavoidable. This will lead to the formation of a silicate-like layer close to the Interface. For a very large range of Hf chemical potential it is very likely that incoming silicon atoms from the Si substrate will form a (HfO2)x(SiO2)y layer, as supported by experimental observations lee . The only way to prevent the formation of a silicate layer would be through the use of an oxygen deficient HfO2, as has been observed Lim . However, even if this silicate formation is prevented in the initial steps of dielectric growth, it seems almost impossible that it will not be formed given the necessary further thermal processing steps, unless a barrier for the Si diffusion being introduced Bastos .

###### Acknowledgements.

This research was supported by Brazilian agencies FAPESP and CNPq. We thank CENAPAD-SP for computer time.


text:
# More about the Laughlin droplet

## Abstract

We compare the exact diagonalization ground state wave function (calculated without any restriction) of a two-dimensional droplet in a perpendicular magnetic field with the Laughlin ansatz, as the number of electrons increases. The fully spin-polarized case for filling factor $`1/3`$ without lateral confining potential or Zeeman effect is considered. We observe that the overlap decreases as the number of electrons increases.

Quantum Hall effect, quantum dot, incompressible (magic) states, exact diagonalization.

Most of the knowledge and understanding of the fractional quantum Hall effect states is due to the existence of Laughlin’s wave function lau , a simple and intuitive analytical expression that matches very well with the ground states (gs) of the system characterized by the filling factor $`\nu =1/q`$, where $`q`$ is an odd integer. This has been tested by numerical calculations of finite systems for different geometries hal ; hal2 ; he ; rez  as well as experimental results related with measurable quantities which can be derived from its properties, as is the case of the fractional charge of some of the excited states gol . We analyze the overlap between this ansatz and the gs wave function that comes out from numerical exact diagonalization of the Hamiltonian for finite systems. We also consider other properties of the gs.

Laughlin’s wave function for $`\nu =1/q`$ is given by:

$$\mathrm{\Psi }_L=\underset{i<j}{}(z_iz_j)^qe^{{\scriptscriptstyle z_i^2/4}}$$
(1)
where $`z=x+iy`$ (in units of $`l_B=\sqrt{\frac{\mathrm{}c}{eB}}`$, $`B`$ being the magnetic field) and where the symmetric gauge has been assumed. The physical meaning is easily deduced from its form: the minimization of the energy in the partially filled lowest Landau level is obtained when each electron sees zeros of order $`q`$ at the positions of the other electrons in such a way that there are non free zeros which would increase the energy. It can be verified that if the polynomial part is expanded, one obtains a combination of Slater determinants built up from Fock Darwin jac  single particle wave functions of the type:

$$\varphi _{sp}z^me^{z^2/4}$$
(2)
(solutions of the non-interacting system) where $`m`$ is the single particle angular momentum. The expansion is a homogeneous polynomial of degree $`M`$ (the total angular momentum) and contains single particle angular momenta up to a maximum given by,

$$m_L=q(N1)$$
(3)
where $`N`$ is the number of electrons. In order to analyze the evolution of the wave function given by Eq.(1) with $`N`$, we consider the cases $`N=2`$, $`3`$ and $`4`$ keeping the condition given by $`\nu =1/3`$ . We compare Eq.(1) with the wave function that comes out from the exact diagonalization (we only considered the fully spin polarized case). In the exact calculation no Zeeman or kinetic energy contributions are considered, electrons are confined within a disk due to the restriction on the total angular momentum.

For $`N=2`$ the dimension of the subspace characterized by the well defined quantum numbers $`M=3`$, the angular momentum needed to obtain $`\nu =1/3`$ for $`N=2`$ ($`M=qN(N1)/2`$) and total spin along the z-direction $`S_z=1`$ is equal to two, that is to say, the exact gs is a combination of two Slater determinants in which the largest possible single particle angular momentum is given by,

$$m_{ex}=M\frac{(N1)(N2)}{2}$$
(4)
or $`m_{ex}=3`$ in this case, in agreement with the value obtained from Eq.(3). Furthermore, the linear combination of Slaters is the same in both cases (in Eq.(1) and in the diagonalization) and consequently the overlap is equal to one. That is to say, for $`N=2`$ the Laughlin wave functiuon is exact. For $`N=3`$ ($`M=9`$) the situation is not the same: the number of Salters involved in the g.s. is not the same, the maximum $`m`$ is also different and as a consequence the overlap is lower than one. This tendency is increased for $`N=4`$. A summary of the results is given in Table 1 and Figs.1 and 2 below. $`N_L`$ and $`N_{ex}`$ are the total number of Slater determinants included in the Laughlin and exact ground states respectively, i.e., $`N_{ex}`$ is the dimension of the Hilbert space characterizad by $`M`$ and $`S_z`$. $`P_L`$ and $`P_{ex}`$ are the normalized weights of the most important Slater within the expansion of the g.s. respectively. This Slater determinant contains, for all the values of $`N`$ considered, a packet of successive single particle angular momenta separated from the center of the dot in such a way that they produce a compact ring. Remarkably this compact structure made up of successive single particle angular momenta wavefunctions ( $`m=2,3`$ and $`4`$ for $`N=3`$) differs from the structure intuitively suggested by that of electrons surrounded by magnetic quantum fluxes. Namely, the zeros at the electronic positions seem to emphasize the short range character of the interaction rather than the separation between electrons. Fig.1 shows the overlap $`S=\mathrm{\Psi }_L\mathrm{\Psi }_{ex}^2`$ for $`N=2`$, $`3`$ and $`4`$ and Fig.2 shows the values of $`P_L`$ and $`P_{ex}`$ for different N. The tendency seems to indicate that the overlap will decrease as $`N`$ increases and so, as the function given by Eq.(1) is a especially good approximation for a low number of electrons, it must be taken with some care for large $`N`$ in a finite system.

It must be emphasized that the tendency of the overlap to worsen as the number of electrons increases would be a trivial result if the gs were a single Slater determinant and also if we used differet single particle wavefunctions: the exact wavefunctions for $`\mathrm{\Psi }_{ex}`$ and approximate wavefunctions for $`\mathrm{\Psi }_L`$. In this case, the overlap would scale as $`(1ϵ)^N`$ $`ϵ`$ being a small number for good trial functions. However this is not our case: the gs is a linear combination of several Slater determinants and in addition, we use the same Slaters to build up the linear combinations in each case, that is to say, the overlap depends on the coefficients, namely, on the electron-electron interaction and is not directly related to $`N`$. As a consequence, for example, Laughlin’s ansatz is the exact result for filling factor $`\nu =1`$ independently of $`N`$, since only one determinant is involved.

Special attention must be paid to the increasing difference between $`N_L`$ and $`N_{ex}`$ due to the fact that the ratio between the weights of the Slater determinants that are lacking in $`\mathrm{\Psi }_L`$ related to the weights of the Slaters which are included is not negliglible. As an example, for $`N=4`$ one of the Slaters lacking in the $`\mathrm{\Psi }_L`$ has a normalized weight of $`0.002`$ compared with that of the most important one which is $`0.331`$. Furthermore, $`\mathrm{\Psi }_L`$ loses the small and the large values of $`m`$ and so the center as well as the edge of the droplet are poorly reproduced. As an example for $`N=5`$, there are $`27`$ Slaters that contain $`m=0`$ and only $`16`$ of them are included in the expansion of $`\mathrm{\Psi }_L`$. However, in spite of the fact that the overlap of the gs wavefunction and Laughlin’s ansatz worsens as $`N`$ increases, it is not the case for some expected values, especially for the energy. If $`H`$ has only the Coulomb contribution (as is appropriate for a fully spin polarized system in the lowest Landau level regime), and if we define the discrepancy $`D`$ as $`[(E_{ex}E_L)/E_{ex}]\times \mathrm{\hspace{0.17em}100}`$ where $`E_{ex}=\mathrm{\Psi }_{ex}H\mathrm{\Psi }_{ex}`$ and $`E_L=\mathrm{\Psi }_LH\mathrm{\Psi }_L`$ the results obtained are:

For $`N=2`$, $`D=0`$, $`D=10.5`$ for $`N=3`$ and $`D=1.3`$ for $`N=4`$, i.e. improving the result as $`N`$ increases.

There are several previous studies on the Laughlin wave function for finite systems which devoted special attention to the study of the edge states: Mitra and MacDonald mit  have analyzed the angular momentum distribution function for a droplet and found that the occupations are peaked at the edge (for $`N=15,20`$ and $`25`$) and that it has a rapid decline. We believe that it can be a consequence of the reduced base implied in $`\mathrm{\Psi }_L`$ as it loses the small and large values of the single particle angular momenta. Tsiper and Goldman tsi  compare the density of a droplet obtained from Eq.(1) and from exact diagonalization for $`N=5`$ to $`12`$ and $`\nu =1/3`$. They obtain important differences at the centre of the droplet and a nearly exact coincidence at the edge. The difference in the electron-electron interaction implied in the Laughlin ansatz (short range interaction) and in the exact calculation (Coulomb interaction) is invoked in order to explain the difference and the formation of striped states. We believe that their result could be a consequence of the procedure used to obtain the results: within the exact calculation they truncate the base of the Hilbert space, i.e., the $`m_{ex}^{^{}}`$ considered is obtained by increasing Hilbert space until overlap $`S`$ converges to at least three significant digits. However this procedure forces a precise coincidence with $`\mathrm{\Psi }_L`$ at the edge, giving no information about the total weight of the rest of the members of the base of the Hilbert space, which is not necessarily negligible, as we mentioned previously. As an example, for $`N=12`$ they consider $`m_{ex}^{^{}}=35`$ ($`m_L=33`$) while $`m_{ex}=143`$ is the exact single-perticle maximum angular momentum involved. At the center, their numerical calculation contains all the Slaters without restriction and they obtain strip-like oscillations on the radial electron densities which has been proved to be responsible for the observed unexpected behavior of the current-voltage power law yan ; wan .

Finally we conclude that for finite systems the overlap between the Laughlin wave function and the exact results is worse as $`N`$ increases. This result is similar to that obtained previously by Haldane hal  for spherical geometry that mimics a two dimensional homogeneous system. His conclusion was that for $`N=3`$ the Laughlin type function is the exact solution ( even for Coulomb interaction) but it is not for $`N4`$. In a recent paper by Yannouleas and Landman yan1  a systematic study of a system of $`6`$ electrons in a range of filling factors from $`\nu =1/5`$ to $`\nu =1/9`$ is reported. They conclude that the analytical model of collectively rotating electron molecules (REM) yan  provides better representation of the system. Other references hal ; he  have tested some results obtained with Eq.(1) for some particular values of $`N`$, however our aim is to study the tendency as $`N`$ increases. The contribution of our report refers to the edge as well as the central properties of finite systems which can be appreciably different from those properties obtained by the use of the Laughlin droplet for large $`N`$ and as a consequence, the differences can be significant at the thermodynamic limit. However, the evolution of the overlap of the wave functions does not necessarily characterize the evolution of the expected values of some operators as was previously pointed out for the eneregy operator.

We gratefully acknowledge F. Salvat and J. Soto for helpful discussions. This work has been performed under Grants No. BFM2002-01868 from DGESIC (Spain) and No. 2001GR-0064 from Generalitat de Catalunya.


text:
# Correlations in Ballistic Processes

## Abstract

We investigate a class of reaction processes in which particles move ballistically and react upon colliding. We show that correlations between velocities of colliding particles play a crucial role in the long time behavior. In the reaction-controlled limit when particles undergo mostly elastic collisions and therefore are always near equilibrium, the correlations are accounted analytically. For ballistic aggregation, for instance, the density decays as $`nt^\xi `$ with $`\xi =2d/(d+3)`$ in the reaction-controlled limit in $`d`$ dimensions, in contrast with well-known mean-field prediction $`\xi =2d/(d+2)`$.

Ballistic-controlled reaction processes CPY ; Jiang ; brl ; mp ; bklr ; bkr ; emm ; F ; PRE ; Ley  exhibit rich atypical behaviors, e.g. the persistent dependence of decay exponents on the spatial dimension $`d`$ implying absence of the upper critical dimension; not surprisingly, ballistic-controlled processes proved very challenging to theoretical treatments. The key such process is ballistic aggregation CPY  in which particles merge upon collisions so that mass and momentum are conserved (energy is necessarily lost). This model arises in various contents, e.g., it mimics the merging of coherent structures (like vortices or thermal plumes) and accumulation of cosmic dust into planetesimals Z . The one-dimensional (1D) version has also an interesting connection with dynamics of shocks representing solutions of the inviscid Burgers equation B ; K . Ballistic aggregation was first investigated in a pioneering paper CPY  by Carnevale, Pomeau, and Young who argued that basic physical quantities behave algebraically in the long time limit, e.g., the density decays as

$$n(t)t^\xi ,\xi =\frac{2d}{d+2}$$
(1)
in $`d`$ dimensions. To understand this result, one can use bklr  a rate equation $`dn/dt=n/\tau `$ for the density. The mean time $`\tau `$ between collisions related to the root mean squared (rms) velocity $`V`$, radius $`R`$, and density through $`nV\tau R^{d1}1`$. Mass conservation implies that the average mass is $`mn^1`$. Therefore $`Rn^{1/d}`$ and

$$\frac{dn}{dt}=n^2VR^{d1}=n^{1+1/d}V.$$
(2)
The particle of mass $`m`$ is formed from $`m`$ original particles (we measure mass in units of the initial mass and velocity in units of the initial rms velocity). Assuming velocities of those original particles uncorrelated we find that the average momentum $`p`$ and velocity $`V`$ scale as

$$pm^{1/2},V=p/mn^{1/2}.$$
(3)
Plugging (3) into (2) and solving for $`n(t)`$ yields (1).

Surprisingly the prediction $`\xi =2d/(d+2)`$ for the decay exponent — perhaps the most known result in the field of ballistic-controlled processes — is erroneous. It turns out that the mean-field assumption that velocities of original particles contained within a typical aggregate particle are uncorrelated is incorrect in any finite dimension — only when $`d\mathrm{}`$ and velocities are orthogonal to each other with probability one, they are indeed uncorrelated. The failure of the mean-field no-correlation assumption (3) has not been appreciated because the resulting formula $`\xi _d=2d/(d+2)`$ is correct both for $`d=1`$ and $`d=\mathrm{}`$. (No trivial explanation of the former assertion is known yet the relation to the Burgers equation via the particles$``$shocks mapping B ; K  and the $`t^{2/3}`$ growth of the separations between adjacent shocks established by Burgers many years ago B  prove that $`\xi _1=2/3`$). Since $`\xi _d`$ monotonously increases with $`d`$, it is not surprising that the actual values are not so different from the mean-field prediction (1). Therefore the observed disagreement in two dimensions emm  could be attributed to insufficient scale of the simulations. Interestingly, the beauty of ballistic aggregation in 1D where the model admits an exact solution mp ; F  and exhibits a deep connection to the Burgers equation has supported the incorrect prediction (1) in higher dimensions.

The purpose of this article is twofold. First, we clarify the role of velocity correlations in the general case, where they lead to significant deviations from mean-field predictions. Second, we propose a procedure that allows an analytical treatment of correlations for virtually any ballistic reaction process in the reaction-controlled limit; in particular, this method gives exact decay exponents.

The no-correlation assumption is generally wrong for all ballistic-controlled processes, so we first demonstrate this assertion for one particularly simple process. We choose a toy ballistic aggregation model in which all particles are identical and when two particles moving with velocities $`𝐯_1`$ and $`𝐯_2`$ collide, they form an aggregate particle moving with velocity $`𝐯=𝐯_1+𝐯_2`$. Compared to the original ballistic aggregation model, the toy model has a number of advantageous properties. First, the volume fraction decays indefinitely thereby driving the system into the dilute limit and justifying ignoring multiple collisions. Second, the mean free path $`n^1`$ grows faster than the inter-particle distance $`n^{1/d}`$ for $`d>1`$. These two features indicate that for $`d>1`$ the Boltzmann equation approach is exact at large times.

For the toy model, (2) becomes $`dn/dt=n^2V`$ and the supposed absence of correlations gives $`Vn^{1/2}`$. Thus the mean-field argument implies $`nt^\xi `$ with $`\xi =2`$ independently on dimension $`d`$. Numerically we find that this universality does not hold: $`\xi `$ increases with dimension and approaches the mean-field prediction only when $`d\mathrm{}`$. For instance, we find (with an accuracy better than 1%)

$$\xi \{\begin{array}{cc}1.33\hfill & \text{when }d=1\text{,}\hfill \\ 1.55\hfill & \text{when }d=2\text{,}\hfill \\ 1.65\hfill & \text{when }d=3\text{.}\hfill \end{array}$$
(4)
These results were obtained by solving the Boltzmann equation describing the toy model

$`{\displaystyle \frac{P(𝐯,t)}{t}}`$ $`=`$ $`{\displaystyle 𝑑𝐮𝑑𝐰P(𝐮,t)P(𝐰,t)|𝐮𝐰|\delta (𝐮+𝐰𝐯)}`$ (5)
$``$ $`2P(𝐯,t){\displaystyle 𝑑𝐰P(𝐰,t)|𝐯𝐰|}.`$
We have solved this equation numerically implementing a Direct Monte Carlo (DMC) simulation scheme (see e.g. Bird  for the general method, and PRE  for an application to a ballistic-controlled reaction process). The idea is to rephrase Eq. (5) as a stochastic process. In each step two particles, say with velocities $`𝐮`$ and $`𝐰`$, are selected at random among a population of $`N`$ particles, and the reaction happens with a probability proportional to $`|𝐮𝐰|`$. If the reaction has been accepted, a new particle of velocity $`𝐮+𝐰`$ replaces two original particles, so the number of particles changes to $`N1`$. The time is incremented by $`(N^2|𝐮𝐰|)^1`$, and the process is iterated again. This numerical scheme allows us to treat systems with initial number of particles of the order of $`10^7`$. The master equation associated to this Markov chain is precisely (5) so that we obtain the numerically exact solution of our problem. The exponent values (4) significantly differ from the mean-field prediction $`\xi =2`$, and leave no doubt that the no-correlation assumption is wrong.

In the Boltzmann equation (5), the relative velocity $`|𝐯𝐰|`$ gives the rate of collisions and its non-linear character makes analytical progress hardly possible. An old trick to overcome this difficulty is to replace the actual relative velocity by the rms velocity max . This results in the Maxwell model that played an important role in the development of kinetic theory classical ; e . For the toy model, we have (hereafter the dependence on time is suppressed for ease of notation)

$$\frac{1}{V}\frac{P(𝐯)}{t}=𝑑𝐮P(𝐮)P(𝐯𝐮)2nP(𝐯).$$
(6)
Integrating (6) we find that the density $`n=𝑑𝐰P(𝐰)`$ satisfies $`dn/dt=n^2V`$ while $`nV^2=𝑑𝐰w^2P(𝐰)`$ remains constant. Hence $`V=n^{1/2}`$ and $`\xi =2`$ showing that the mean-field no-correlation approach is essentially the Maxwell model in context of ballistic processes bk . The Maxwell model is an uncontrolled approximation to the Boltzmann equation for the hard sphere gas and, not surprisingly, the exponents found within this approach are generally erroneous (see Rqvhp  for an alternative simplification, the so-called very hard particle approach). Of course, one could anticipate that the exponent $`\xi =2`$ characterizes the Maxwell model without computations — the essence of the Maxwell model, that is the fact that collisions are completely random, assures that the no-correlation condition does hold.

We now present an argument that emphasizes the role and importance of correlations between velocities of colliding particles and applies to all ballistic-controlled reaction processes. The key point is to supplement an evolution equation for the mass density by an evolution equation for the density of kinetic energy. For an arbitrary ballistic-controlled reaction process we denote $`P(m,𝐯,t)`$ the joint mass-velocity distribution function, and $`e=mv^2`$ the kinetic energy of a given particle (for the toy model, we set $`m=1`$). The evolution equations for the density $`n`$ and kinetic energy density $`nE=mv^2P(m,𝐯,t)𝑑m𝑑𝐯=nmv^2`$ read

$`{\displaystyle \frac{dn}{dt}}={\displaystyle \frac{n}{\tau }},{\displaystyle \frac{d(nE)}{dt}}={\displaystyle \frac{n\mathrm{\Delta }e_{\text{coll}}}{\tau }}.`$ (7)
The first equation is just the definition of the time dependent collision frequency $`1/\tau `$, the second additionally contains the kinetic energy $`\mathrm{\Delta }e_{\text{coll}}`$ lost on average in a binary collision. In the scaling regime the quantities $`\mathrm{\Delta }e_{\text{coll}}`$ and $`E=mv^2`$ exhibit the same time dependence, so the dissipation parameter $`\alpha =\mathrm{\Delta }e_{\text{coll}}/E`$ is asymptotically time independent. From Eqs. (7) we get $`d\mathrm{ln}(nE)/d\mathrm{ln}n=\alpha `$, or $`V^2=nEn^\alpha `$. The mean free path argument $`\tau ^1nVR^{d1}t^1`$ gives $`n^{1/d}Vt^1`$ for ballistic aggregation. Combining these two relations and the definition of $`\xi `$ we obtain $`\xi =(1/d+\alpha /2)^1`$. Similarly for the toy model $`nV^2n^\alpha `$ and $`nVt^1`$ leading to $`\xi =2/(1+\alpha )`$.

To use this formalism, we must precisely define the collisional average involved in (7). An average change of a quantity $`𝒜(1,2)`$ in a binary collision is radius

$$\mathrm{\Delta }𝒜_{\text{coll}}=\frac{d1d2|𝐯_1𝐯_2|^\nu [\mathrm{\Delta }𝒜(1,2)]P(1)P(2)}{d1d2|𝐯_1𝐯_2|^\nu P(1)P(2)},$$
(8)
where we have used shorthand notations $`i=(m_i,𝐯_i)`$ and $`di=dm_id𝐯_i`$ ($`i=1,2`$). In the key case of hard spheres we have $`\nu =1`$, whereas the cases $`\nu =0`$ and $`\nu =2`$ correspond the Maxwell and very hard particle models Rqvhp , respectively. We now illustrate the formalism for the toy model. The kinetic energy lost in a collision is $`\mathrm{\Delta }e=𝐯_1^2+𝐯_2^2(𝐯_1+𝐯_2)^2=2𝐯_1𝐯_2`$. Hence

$`\mathrm{\Delta }e_{\text{coll}}=2{\displaystyle \frac{𝑑𝐯_1𝑑𝐯_2|𝐯_1𝐯_2|^\nu (𝐯_1𝐯_2)P(𝐯_1)P(𝐯_2)}{𝑑𝐯_1𝑑𝐯_2|𝐯_1𝐯_2|^\nu P(𝐯_1)P(𝐯_2)}}.`$
For the Maxwell model ($`\nu =0`$), the isotropy of $`P(𝐯,t)`$ shows that $`\mathrm{\Delta }e_{\text{coll}}=0`$, so $`\alpha =0`$ and $`\xi =2/(1+\alpha )=2`$ in agreement with our previous calculation. Similarly for very hard particles ($`\nu =2`$, see Rqvhp ) we use isotropy to simplify $`\alpha `$ and arrive at

$`\alpha =2{\displaystyle \frac{𝑑𝐯_1𝑑𝐯_2(𝐯_1𝐯_2)^2P(𝐯_1)P(𝐯_2)}{\left[𝑑𝐯𝐯^2P(𝐯)\right]^2}}.`$
The isotropy allows to compute the ratio of the integrals to yield $`\alpha =2/d`$ leading to $`\xi =2d/(d+2)`$. For other values of $`\nu `$, including the case of interest $`\nu =1`$, the dissipation parameter $`\alpha `$ depends on details of the velocity distribution, and isotropy alone is not sufficient to determine $`\alpha `$. The reason for the failure of the mean-field argument — which amounts to the complete neglect of collisional correlations ($`𝐯_1𝐯_2_{\text{coll}}=0`$) — is now clear: in general, a collision involving a pair $`(𝐯_1,𝐯_2)`$ with a negative product $`𝐯_1𝐯_2<0`$ has a higher probability than a collision with $`𝐯_1𝐯_2>0`$. The dissipation parameter $`\alpha `$ is therefore positive so that $`\xi =2/(1+\alpha )<2`$. Thus the mean-field prediction $`\xi =2`$ is an upper bound for $`\xi `$.

The above framework applies to any irreversible process with ballistic transport. For ballistic aggregation, the omission of collisional correlations amounts to setting $`\alpha =1`$, i.e. that the typical energy dissipated in a collision is the mean kinetic energy per particle. However, particles with larger velocities undergo more frequent collisions so that the mean energy dissipated exceeds the energy of a typical particle. Hence $`\alpha =\mathrm{\Delta }e_{\text{coll}}/E>1`$ so that $`\xi `$ is smaller than the mean-field prediction $`2d/(d+2)`$. Previous Molecular Dynamics (MD) simulations have shown that $`\xi 0.85\pm 0.04`$ in 2D for low volume fractions, with scaling laws extending over 2 decades in time emm . The DMC technique allows to reach much larger time scales. Figure 1 shows that after an initial transient, the density exhibits a clear power law behavior over 5 decades in time. We estimate $`\xi 0.86\pm 0.005`$, in agreement with MD simulations. The inset displays the behavior of $`E=mv^2`$, the quantity that is (asymptotically) time independent according to the mean-field prediction (3); we find $`Et^{0.28}`$ R . We have also performed DMC and MD simulations in 3D giving $`\xi 1.06\pm 0.01`$. As expected, the actual values of $`\xi `$ are smaller than the mean-field prediction $`\xi =2d/(d+2)`$note1 .

Ballistic-controlled processes are generally intractable analytically. Following the fruitful line of attack on difficult problems — generalize them! — let us consider a process in which colliding particles react with probability $`ϵ`$ and scatter elastically with complementary probability $`1ϵ`$. The mean-field no-correlation argument is so general that it applies to these processes; in particular, according to mean-field the exponent $`\xi `$ is independent on $`ϵ`$. Remarkably, we can now compute the exponent $`\xi `$ for one special value of $`ϵ`$, viz. for $`ϵ0^+`$. In this reaction-controlled limit particles undergo mostly elastic collisions. Therefore, the particles are always at equilibrium, i.e., the velocity distribution is Maxwellian. This key feature makes the problem tractable. Consider for instance the toy model. One can compute

$`\mathrm{\Delta }e_{\text{coll}}=2{\displaystyle \frac{𝑑𝐯_1𝑑𝐯_2|𝐯_1𝐯_2|^\nu (𝐯_1𝐯_2)P(𝐯_1)P(𝐯_2)}{𝑑𝐯_1𝑑𝐯_2|𝐯_1𝐯_2|^\nu P(𝐯_1)P(𝐯_2)}}`$
for arbitrary $`\nu `$ when $`P(𝐯)`$ is Maxwellian help . In particular for the important case $`\nu =1`$ we obtain $`\alpha =1/d`$, so that $`\xi =2d/(d+1)`$. This exact result provides a useful check of numerical scheme (our DMC simulations in two and three dimensions are indeed in excellent agreement with the theoretical prediction). In contrast, the mean-field no-correlation argument predicts $`\xi =2`$ irrespective of the value of $`ϵ`$. We see again that this is correct only in the $`d\mathrm{}`$ limit. Note also that $`\xi =2d/(d+1)`$ which is exact for the reaction-controlled ($`ϵ0^+`$) version of the toy model provides a better ‘guess’ for $`\xi `$ in the original ($`ϵ=1`$) model than the mean-field approach \[compare $`\xi =1`$, 4/3 and 3/2 in 1D, 2D, and 3D to the numerical values (4)\].

Remarkably, the exponent $`\xi `$ in the reaction-controlled limit of ballistic aggregation can be computed even though the mass distribution $`\mathrm{\Pi }(m)=n^1𝑑𝐯P(m,𝐯)`$ is unknown. The important point is that when $`ϵ0^+`$, the joint mass/kinetic energy distribution function factorizes. Then one finds $`\alpha =1+1/d`$, or equivalently $`\xi =2d/(d+3)`$ independently on $`\mathrm{\Pi }(m)`$ help2 . This exact result of course agrees with DMC simulations. Interestingly, it also provides a reasonable approximation of $`\xi `$ for the original ($`ϵ=1`$) aggregation model: $`\xi =0.8`$ in 2D and 1 in 3D, to be compared to 0.86 and 1.06, respectively.

Many other ballistic-reaction processes are solvable in the reaction-controlled limit. For instance for ballistic annihilation brl , there is no exact solution in any dimension yet in reaction-controlled limit, the exact value of the density decay exponent is given by $`\xi =4d/(4d+1)`$. This result is in surprisingly good agreement with numerical values for $`ϵ=1`$: $`\xi =4/5`$ vs. 0.804 Ely  in 1D; $`\xi =8/90.89`$ vs. 0.87 PRE  in 2D; $`\xi =12/130.92`$ against 0.91 PRE  in 3D. We have studied several other ballistic-reaction processes prep , e.g., a simplified ballistic aggregation model in which mass and momentum are conserved yet the radius does not grow. For this model, the mean-field prediction is $`\xi =2/3`$ independently on dimension $`d`$, whereas in the reaction-controlled limit, we get the exact result $`\xi =2d/(3d+1)`$ (i.e. 0.571 in 2D and 0.6 in 3D). It is again instructive to compare these values with numerical results for $`ϵ=1`$: $`\xi 0.60`$ in 2D, and $`\xi 0.62`$ in 3D.

We have shown that correlations between velocities of colliding particles govern the behavior of all reacting processes with ballistic transport. We illustrated importance of correlations on several models and demonstrated that ignoring correlations is equivalent to using the Maxwell model, which is an uncontrolled approximation of the hard-sphere gas. We also devised a procedure that clarifies the role of correlations in the general case and allows an exact computation of decay exponents in the reaction-controlled limit, when particles undergo mostly elastic collisions and therefore are always near equilibrium. The failure of mean-field theory to describe this limit emphasizes inevitable presence of correlations in all reacting processes with ballistic transport.


text:
# Big Bang Nucleosynthesis and the Missing Hydrogen Mass in the UniverseTalk given at CIPANP2003. To be published by the American Institute of Physics in the Proceedings of the Conference on the Intersections of Particle and Nuclear Physics (CIPANP2003) held at Grand Hyatt Hotel – New York City, May 19-24, 2003.

## Abstract

It is proposed that when the era of the big-bang nucleosynthesis ended, almost all of the 75 percent of the observed total baryonic matter remained in the form of hydrogen and continued to exist in the form of protons and electrons. They are present today as baryonic dark matter in the form of intergalactic hydrogen plasma. To test our hypothesis we have investigated the effects of Thomson scattering by free electrons on the reported dimming of Type Ia supernovae. The quantitative results of our calculation suggest that the dimming of these supernovae, which are dimmer than expected and hence more distant than predicted by Hubble expansion, is a result of Thomson scattering without cosmic acceleration.

Recent observations \[1-2\] of Type Ia supernovae (SNe Ia) appear to suggest that the universe is accelerating. The basic idea proposed to account for the acceleration is dark energy or quintessence \[3-6\]. However, as yet there is no direct confirmation of their existence or exact nature. The present work examines whether the question of acceleration can be resolved within the limits of established laws of physics. Consequently, we have investigated the effects of Thomson scattering of photons by free electrons present in the form of H-plasma while propagating from the supernova source to the point of observation. This particular scattering process has been chosen for the following reasons: the scattering cross-section is pure elastic in nonrelativistic region, the total cross-section is a universal constant and it is independent of the incident frequency . Hence characteristics of the atomic spectra, which are relevant in the present work, remain unchanged. The calculation is performed within the framework of Friedmann-Robertson-Walker (FRW) cosmology  for the special case of a flat universe, consistent with the recent cosmic microwave background anisotropy measurements  indicating a spatially flat, critical density universe with $`\mathrm{\Omega }=1`$ and with the inflationary model \[11-13\] of cosmology. Hence in the present investigation we consider a universe, for large-scale structure, consisting only of matter and negligible radiation without cosmological constant in accord with all popular cosmological models prior to the late 1990s. Therefore in terms of standard notation,

$$\mathrm{\Omega }=\frac{\rho (z)}{\rho _c(z)}=1$$
(1)
and

$$\rho _c(z)=\frac{3H^2(z)}{8\pi G}=\frac{3H_0^2}{8\pi G}(1+z)^3$$
(2)
where $`\rho `$ is the total mass density including radiation, $`\rho _c`$ is the critical density, $`G`$ is the gravitational constant, $`H`$ is the Hubble constant and in which we have incorporated the flat universe relation between $`H(z)`$ and the current Hubble constant $`H_0`$.

In light of the above considerations, the total matter density of universe consists of ordinary baryonic matter, baryonic dark matter and nonbaryonic dark matter; thus we have $`\mathrm{\Omega }=\mathrm{\Omega }_{bm}+\mathrm{\Omega }_{bdm}+\mathrm{\Omega }_{nbdm}`$ where each term is in units of $`\rho _c`$, and where the subscripts refer, respectively, to the three aforementioned components. The observed mass density $`\mathrm{\Omega }_m0.30\pm 0.10`$ as estimated from its gravitational pull on visible matter and is assumed to consist of ordinary baryonic matter and baryonic dark matter instead of exotic dark matter ; thus $`\mathrm{\Omega }_m=\mathrm{\Omega }_{bm}+\mathrm{\Omega }_{bdm}`$. Big bang nucleosynthesis and synthesis of heavy elements in stars have been extensively investigated for more than half a century. The quantitative results of atomic abundances of various groups of elements \[15-16\] by mass-fraction of the total are hydrogen $`0.75`$ and rest of the elements $`0.25`$. From the above value of $`\mathrm{\Omega }_m`$ we obtain the magnitude of $`\mathrm{\Omega }_{bm}`$ which includes all the elements formed during big bang nucleosynthesis, and of $`\mathrm{\Omega }_{bdm}`$ which remained as free hydrogen (in the form of mostly protons and equal number of electrons, neutrons having decayed into protons and electrons) within half an hour after the era of big bang nucleosynthesis ended. These magnitudes are $`0.05\mathrm{\Omega }_{bm}0.10`$ and $`0.15\mathrm{\Omega }_{bdm}0.30`$. From this analysis we conclude that when the era of the big bang nucleosynthesis ended, the universe continued expanding for several thousand years or so until the temperature dropped low enough to form neutral atoms. We believe that it is in this time interval that most of the free protons and electrons (contained in $`\mathrm{\Omega }_{bdm}`$) escaped into cosmic space and are most likely present in the form of an intergalactic hydrogen plasma . They are dark because they cannot emit light.

The intensity of the radiation lost to Thomson scattering depends critically on the density of free electrons in the path from the source to the observer. We take the free electron number density at redshift $`z`$ as

$$n(z)=\frac{\mathrm{\Omega }_{bdm}\rho _c(z)}{m_h}$$
(3)
where $`\mathrm{\Omega }_{bdm}\rho _c`$ is the hydrogen mass density and $`m_h`$ is the hydrogen mass. We consider radiation emitted by a Type Ia supernova at redshift $`z_S`$ and received by an infinitesimal volume element of length $`cdt`$ at redshift $`z`$. The fractional reduction of intensity owing to Thomson scattering by free electrons within this volume is

$$\frac{dI}{I(z)}=\sigma _Tn(z)cdt$$
(4)
where $`I(z)`$ is the incident intensity and $`\sigma _T`$ is the total Thomson scattering cross-section by a free electron. The total attenuation by scattering is the integral of eq. (7) over the entire path taken by the light from the supernova. We perform the integration by identifying $`dt`$ with $`dT_H(z)`$ where $`T_H(z)`$ is the Hubble time $`H^1(z)`$. Thus

$$dt=dT_H(z)=\frac{3}{2}H_0^1(1+z)^{5/2}dz.$$
(5)
Combining eqs. (2)-(5) then yields

$$\frac{dI}{I}=\sigma _T\mathrm{\Omega }_{bdm}\frac{9cH_0}{16\pi Gm_h}(1+z)^{1/2}dz.$$
(6)
Integrating over a path from the source to the observer at $`z=0`$ yields

$$I(0)=I_S\mathrm{exp}\left\{\sigma _T\mathrm{\Omega }_{bdm}\frac{3cH_0}{8\pi Gm_h}[(1+z)^{3/2}1]\right\}$$
(7)
where $`I_S`$ is the light intensity at the source and in which the exponential factor represents the loss of intensity by Thomson scattering.

Observations of redshift and distance modulus $`mM`$ for Type Ia supernovae reported by Riess et al.  are listed in Table 1 in columns 1 and 2. Column 3 contains the distance computed for each distance modulus according to

$$R_{\text{obs}}=10^{(mM25)/5}$$
(8)
where $`R_{\text{obs}}`$ is the distance in Mpc. Corrections for Thomson scattering are effected by multiplying the observed values of $`mM`$ by the ratio $`I(0)/I_S`$ from eq. (7) and columns 4-6 contain the corrected distances as computed from eq. (8); these three columns correspond to the extreme values of the hydrogen mass fraction $`\mathrm{\Omega }_{bdm}`$ listed above and to their midpoint, $`\mathrm{\Omega }_{bdm}0.23`$. The calculations employ $`H_0=65`$ km-s-1-Mpc-1. The last column lists values of the theoretical luminosity distance $`D_L`$ which, for a flat universe with deceleration parameter $`q_0=1/2`$, is given by

$$D_L=\frac{2c}{H_0}(1+z\sqrt{1+z}).$$
(9)
*Conclusion.* The effects of Thomson scattering on enhanced dimming of SNe Ia presented in Table 1 suggests that (i) the recently observed supernovae data can be understood without dark energy; (ii) the amounts of baryonic ordinary matter, baryonic dark matter and nonbaryonic dark matter are in the ranges of 5-10%, 15-30% and 60-80%, respectively; and (iii) the total matter density $`\mathrm{\Omega }`$ consists of 5-10% baryonic ordinary matter and 90-95% total dark matter, consistent with the understanding of most cosmologists prior to late 1990s. Further details of the present investigation will be published elsewhere.

REFERENCES

1. Riess, A.G. et al, *Astron. J.* 116, 1009-1038 (1998).

2. Perlmutter, et al, *Astrophys. J.* 517, 565-586 (1999).

3. Wang, L. and Steinhardt, P. J., *Astrophys. J.* 508, 483-490 (1998).

4. Wang, L., Caldwell, R. R., Ostriker, J. P. and Steinhardt, P. J., *Astrophys. J.*

530, 17-35 (2000).

5. Perlmutter, S., Turner, M. S. and White, M., *Phys. Rev. Lett.* 83, 670-673

(1999).

6. Peebles, P. J. E. and Ratra, B., *Rev. Mod. Phys.* 75, 559-606 (2003).

7. Heitler, W., *The Quantum Theory of Radiation*, 3rd ed.,(Dover, New York,

1984), pp. 34-35.

8. Kolb, E. W. and Turner, M. S., *The Early Universe*, paperback ed., (Addison-

Wesley, Massachusetts, 1994), pp. 47-86.

9. de Bernardis, P. et al., *Nature* 404, 955-959 (2000).

10. Netterfield, C. B. et al., *Astrophys. J.* 571, 604-614 (2002).

11. Guth, A. H., *Phys. Rev.* D23, 347-356 (1981).

12. Linde, A. D., *Phys. Lett.* B108, 389-393 (1982).

13. Albrecht, A. and Steinhardt, P. J., *Phys. Rev. Lett.* 48, 1220-1223 (1982).

14. Bahcall, N. A., Ostriker, J. P., Perlmutter, S. and Steinhardt, P. J.,

*Science* 284, 1-16 (1999).

15. Burbidge, E. M., Burbidge, G. B., Fowler, W. A. and Hoyle, F., *Rev. Mod.*

*Phys.* 29, 547-650 (1957).

16. Boesgaard, A. M. and Steigman, G., *Ann. Rev. Astron. Astrophys.* 23,

319-378 (1985).

17. Weinberg, S., *Gravitation and Cosmology* (Wiley, New York, 1972), p. 500.

18. Peebles, P. J. E., *Principles of Physical Cosmology*, (Princeton Univ. Press,

Princeton, 1993), p.21.