text:
# Bistability and macroscopic quantum coherence in a BEC of ⁷𝐿𝑖.

## I Introduction

The availability of Bose-Einstein condensates of trapped cold atoms RbRef ; Ketterle ; LiRef ; Li97  has opened the possibility of a laboratory engineering of quantum states with a large number of atoms Walsworth  (around a thousand for $`{}_{}{}^{7}Li`$ Li97 ; RbCond ).

One of the most challenging endeavours of quantum engineering is the evidence of superposition states \[so called SC=Schrödinger cat (SC)\] whose mutual interference is called macroscopic quantum coherence (MQC). SC have been observed, e.g., for states of a trapped ion Monroe  and of a microwave field in a high-Q cavity Brune .

In this paper we demonstrate the reliable preparation of SC consisting of a Bose-Einstein condensate (BEC) of $`{}_{}{}^{7}Li`$ atoms that have negative scattering length, trapped in a double-well potential. Realistic calculations have been offered for macroscopic quantum tunneling (MQT) Leg ; Sala . Indeed, combining the kinetic and potential terms of a harmonic trap with the inter-particle attraction yields a metastable state for $`N<N_c`$ ($`N_c=`$ critical population for an attractive BEC). Quantum tunneling from this metastable state towards the collapsed state, which would otherwise be reached for $`N>N_c`$, has been shown to be feasible.

A BEC of atoms with negative scattering length, trapped in a double well, undergoes a space symmetry breaking beyond a threshold number of atoms $`N_i`$, whereby two stable states are formed. This phenomenon has been dealt with theoretically by a two mode approach Ho .

We have studied the problem by finding numerically the stationary solutions of the Gross-Pitaevskii (GP) equation discretized over a space lattice, with reference to the $`{}_{}{}^{7}Li`$ case artipre . Once the stationary solutions have been found, we introduce a quantum two mode model, with the two modes chosen in such a way as to reproduce the stationary solutions of GP. The model shows the feasibility of macroscopic quantum coherence (MQC).

Bistability occurs only for an attractive interatomic potential (negative scattering length). Two proposals for reaching SC in a BEC with repulsive atoms have been put foward Cirac ; Gordon . They both require a Raman coupling between two different components; Both the papers consider copropagating light beams. In Ref. Cirac , MQC is shown to require values of the scattering length between atoms of different magnetic number substantially larger than the scattering length between atoms of equal magnetic number. This requirement is too strong, since no experimental technique is today accessible to provide such a difference; furthermore, if such a difference could be achieved, outstanding symmetry-breaking effects would occur Savage . Reference Gordon  introduces a time-dependent evolution, so that SC is reached over a time of the order of 1 sec. However, no clear-cut experimental test is offered to discriminate between SC and a statistical mixture of two separate states.

## II Bistability and symmetry breaking

We refer to a condensate of $`{}_{}{}^{7}Li`$ atoms trapped in a double-well potential. A suitable model for it is given by

$$V(\stackrel{}{x})=\frac{1}{2}m\left[\omega _{}^2x_1^2+\omega _{}^2(x_2^2+x_3^2)\right]+Acos\left(2\pi \frac{x_1}{\sigma }\right)$$
(1)
The quadratic part is due to the interaction of the atoms with the magnetic field of the trap. According to laboratory implementations Li97  we choose $`\omega _{}=2\pi \times 130s^1`$ and $`\omega _{}=2\pi \times 150s^1`$. The additional term is generated by two opposite laser beams in a standing wave configuration. A suitable choice of the standing-wave parameters yields a double well-potential. Taking into account the interatomic interaction, the atomic system is described by a macroscopic wave function $`\psi `$ that satisfies the GP equation

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}^2\psi +V\psi +g|\psi |^2\psi (H+g|\psi |^2)\psi ,$$
(2)
with $`H=\mathrm{}^2/(2m)^2+V`$. Here, $`m7a.u.`$ is the mass of the lithium atom, and $`g=(4\pi \mathrm{}^2/m)a_s`$, where $`a_s`$ is the s-wave scattering length for $`{}_{}{}^{7}Li`$, $`a_s1.45nm`$. For a small number of atoms the GP nonlinearity can be neglected and Eq. (2) reduces to an ordinary Schrödinger equation. In such a case and for a sufficiently high barrier, the lowest energy level is described by a two-peak wave function symmetric with respect to inversion of the space axes, once the coordinate origin coincides with the trap center (Fig. 1, dashed-dot line).

Figure 1 reports the spatial distribution of the ground state of Eq. (2) for different numbers of trapped atoms; we have used a numerical method that consists in solving GP on a discrete space lattice and evaluating the lowest energy state. The barrier is specified by the two parameters $`A_nA/\mathrm{}=2650s^1`$ and $`\sigma =5\mu m`$ (see Eq. (1)). As shown in the figure, for $`N=450`$ the distribution is symmetric; instead for $`N=500`$ the nonlinear term is sufficient to destabilize the symmetric state, giving rise to two asymmetric stable states. For $`N=600`$ one well is almost empty. As we increase the number of atoms, the nonlinearity plays a relevant role. By a self-consistent argument we realize that the symmetric wave function becomes unstable and we can have two new minimal energy states with distribution no longer symmetric for inversion (symmetry breaking). Indeed, let us assume a distribution as in Fig. 1 (dashed or solid line); then the effective potential for such a distribution, due to the sum of the external potential with $`g|\psi |^2`$ is an asymmetric double well with the lower minimum corresponding to the higher population peak. For a sufficiently high nonlinear term the potential imbalance stabilizes the asymmetric distribution as in Fig. 1.

We confirm the numerical calculation by the following analytic model. Let $`\psi _a`$ be the equilibrium symmetrical wave function (either stable or unstable), and $`\psi _b`$ be a suitable antisymmetrical wave function such that the weighted sum of the two wave functions lowers either one of the two peaks. In the two-dimensional space of these wave functions, any other one can be expressed as

$$\psi (\stackrel{}{x})=a\psi _a(\stackrel{}{x})+b\psi _b(\stackrel{}{x}).$$
(3)
Without loss of generality we can choose $`\psi _a`$ and $`\psi _b`$ as real functions and thus consider real values for $`a`$ and $`b`$. As we take $`\psi _a^2d^3x=\psi _b^2d^3x=1`$, it follows that $`a^2+b^2=\psi ^2d^3xM(a,b)`$, where $`M(a,b)`$, the total number of atoms in the condensate, depends upon $`a`$ and $`b`$. The energy is then given by

$$=a^2H_{aa}+b^2H_{bb}+\frac{1}{2}(a^4I_{aa}+b^4I_{bb})+3a^2b^2I_{ab}$$
(4)
where

$$H_{aa}=\psi _a^{}H\psi _ad^3x,H_{bb}=\psi _b^{}H\psi _bd^3x,$$

$$I_{aa}=g\psi _a^4d^3xx,I_{bb}=g\psi _b^4d^3x,$$

$$I_{ab}=g\psi _a^2\psi _b^2d^3x.$$

We look for the minimal of energy with the constraint of a fixed number $`N`$ of condensed atoms. These constrained minima satisfy the relations

$$\frac{(\mu M)}{a}=\frac{(\mu M)}{b}=0,$$
(5)
where $`\mu `$ is a Lagrange multiplier. We thus solve for $`a`$ and $`b`$, with the condition $`N=a^2+b^2`$. The solutions $`a=0`$ and $`b=0`$ correspond, respectively, to the antisymmetric function and the symmetric one. The other solutions, for both $`a`$ and $`b`$ nonzero, yield $`a^2`$, $`b^2`$ values as functions of $`\mu `$. Using the constraint of fixed $`N`$, we eliminate $`\mu `$ and find

$`a^2={\displaystyle \frac{H_{aa}H_{bb}+N(3I_{ab}I_{bb})}{6I_{ab}I_{aa}I_{bb}}}`$ (6)
$`b^2={\displaystyle \frac{H_{bb}H_{aa}+N(3I_{ab}I_{aa})}{6I_{ab}I_{aa}I_{bb}}}`$ (7)
Since $`g`$ is negative, the denominator is always negative. Indeed, $`\psi _a^2`$ and $`\psi _b^2`$ are almost equal at each point $`x`$. For low $`N`$, the dominant terms in the numerators have opposite sign, thus one of the two squares has to be negative, which means that there is no solution with asymmetric wave function. On the other hand, for sufficiently large $`N`$, the second term in the numerator of the two equations can compensate for the positive one, since the quantities $`3I_{ab}I_{aa}`$, $`3I_{ab}I_{bb}`$ are always negative. Thus, we have proved that two asymmetrical steady states exist beyond a threshold value of $`N`$.

In order to prove that the two states are stable, it is sufficient to show that the symmetrical state becomes unstable above threshold, that is,

$$\frac{d^2}{db^2}|_{b=0}0$$
(8)
where we have $`a^2=Nb^2`$. It is easily found that

$$\frac{d^2}{db^2}|_{b=0}=2(H_{bb}H_{aa})+2N(3I_{ab}I_{aa})$$
(9)
Going back to Eq. (7), the two stationary asymmetrical states occur when the numerator changes sign. In fact, the numerator is the right side of Eq. (9), hence, at the critical point the symmetrical state becomes unstable.

One can easily evaluate the threshold value $`N_i`$ of $`N`$ for which this symmetry breaking occurs.

## III Macroscopic Quantum Coherence

In such a bistable situation the energy displays two minima of equal value in the infinite-dimensional phase space. From a classical point of view, as the system is in its lowest energy state, the condensate is localized in either one of the two minima, where it will remain in the absence of thermal noise once we keep the atom number constant. Since, however, the condensate is a mesoscopic system, quantum fluctuations play a relevant role. This can be shown by second quantization of the field, replacing the c-number macroscopic wave function by field operators.

Quantum fluctuations allow the passage from one to the other stable state without thermal activation, by pure quantum tunneling. Furthermore, due to the coherent nature of the process, we expect coherent oscillations between the two wells, that is, MQC. We now evaluate the tunneling rate as a function of the system parameters showing the feasibility of MQC.

The most natural way of evaluating the tunneling rate consists in finding the two lowest eigenvalues of the Hamiltonian and taking their difference. Indeed, the sum and difference of the corresponding states are respectively the alive and dead states of SC, and the transition time is half the period corresponding to the frequency difference.

The problem is simplified by reducing it to two degrees of freedom by the expansion of Eq.(3).

We select the basis functions $`\psi _a`$ and $`\psi _b`$ as follows. Calling $`\psi _0(\stackrel{}{x})`$ the minimal-energy wave function of the GP problem (see e.g. Fig. 1), we take for $`\psi _a`$ and $`\psi _b`$, respectively, the symmetric and antisymmetric sums $`\psi _0(\stackrel{}{x})\pm \psi _0(\stackrel{}{x})`$. Expansion (3) with these $`\psi _a`$ and $`\psi _b`$ includes the original functions $`\psi _0(\stackrel{}{x})`$ and $`\psi _0(\stackrel{}{x})`$ for suitable values of $`a`$ and $`b`$. Furthermore, it simplifies the form of the Hamiltonian, as we see right now nota .

In second quantization, $`a`$ and $`b`$ in Eq. (3) become the operators $`\widehat{a}`$ and $`\widehat{b}`$, obeying Bose commutation rules with their conjugates $`\widehat{a}^{}`$, $`\widehat{b}^{}`$. Exploiting the operator version of Eq. (3) and its adjoint, the Hamiltonian becomes

$$=\widehat{a}^{}\widehat{a}H_{aa}+\widehat{b}^{}\widehat{b}H_{bb}+\frac{1}{2}(\widehat{a}^2\widehat{a}^2I_{aa}+\widehat{b}^2\widehat{b}^2I_{bb})+$$
(10)
$$\left(\frac{1}{2}\widehat{a}^2\widehat{b}^2+\frac{1}{2}\widehat{b}^2\widehat{a}^2+2\widehat{a}^{}\widehat{b}^{}\widehat{a}\widehat{b}\right)I_{ab}$$

where the coefficients $`H_{aa}`$, $`H_{bb}`$, $`I_{aa}`$, $`I_{bb}`$ and $`I_{ab}`$ are the same as in Eq. (4).

We consider the basis of eigenvectors of the number operators

$$|0,N>,|1,N1>,\mathrm{}|N,0>$$
(11)
where $`\widehat{a}^{}\widehat{a}|k,m>=k|k,m>`$ and $`\widehat{b}^{}\widehat{b}|k,m>=m|k,m>`$. Let us call

$$H_{l,k}=<l,Nl||k,Nk>$$
(12)
the generic matrix element of the Hamiltonian on the above basis. We evaluate the eigenvalues of this matrix. We have considered two different wavelengths and two amplitude values of the applied field.

In Fig. 2 we report on the first two excited energy levels versus the number of condensed atoms beyond the threshold value. In such figure we keep constant the barrier height and just vary the barrier width. Notice that for increasing $`\sigma `$ the maximum tunneling frequency reduces, but the slope at which it reduces for increasing $`NN_i`$ is less steep.

In Fig. 3 we keep fixed the barrier width and change its height. Here too the maximum frequency decreases for increased heights, but again as in Fig. 2, the slope decreases for increasing $`NN_i`$.

In Fig. 4 we reduce the threshold for BEC breaking by increasing $`\omega _{}`$. Precisely, we report on the first two excited states for $`\omega _{}=2\pi \times 600s^1`$, $`A_n=2000s^1`$ and $`\sigma =5\mu m`$. The threshold value is $`N_i190`$.

In the number representation there is no explicit evidence of a SC as a two-peak distribution (see Fig. 5). We look for a suitable observable, whose probability distribution provides such an evidence. We take the first component of the barycenter coordinate $`x_c=(1/N)x_1|\psi |^2d^3x`$ as the appropriate variable since the corresponding classical states \[minima of the Hamiltonian (4)\] have separated barycenters. It is associated with the operator

$$\widehat{x}_c=\frac{1}{N}x_1\psi _a(\stackrel{}{x})\psi _b(\stackrel{}{x})d^3x(\widehat{a}^{}\widehat{b}+\widehat{a}\widehat{b}^{})$$
(13)
Besides a c-number factor, the observable is thus $`\widehat{M}=\widehat{a}^{}\widehat{b}+\widehat{a}\widehat{b}^{}`$. The associated probability density is $`P(m)=|<m|\varphi _0>|^2`$, where $`|\varphi _0>`$ is the ground state of Hamiltonian (10) and $`|m>`$ is the eigenstate of $`\widehat{M}`$ with eigenvalue $`m`$. Notice that once the number of atoms is fixed the eigenvectors are not degenerate.

Since we know the components of $`|\varphi _0>`$ on the number basis we must express $`|m>`$ with respect to that basis, that is,

$$|m>=\underset{k}{}c_k^m|k,Nk>.$$
(14)
Applying the operator $`\widehat{M}`$ to the above ket we have

$$m|m>=\underset{k}{}c_k^m(\widehat{a}^{}\widehat{b}+\widehat{a}\widehat{b}^{})|k,Nk>.$$
(15)
Projecting on $`<l,Nl|`$ the above ket it follows that

$$\underset{k}{}M_{l,k}c_k^m=mc_l^m,$$
(16)
where

$$M_{l,k}=<l,Nl|\widehat{M}|k,Nk>.$$
(17)
Since $`M_{l,k}`$ are known, solving the eigenvector equation (16) we can evaluate the coefficients $`c_k^m`$ and hence $`P(m)`$.

The two-peak distribution $`P(m)`$ is plotted in Fig. 6(a) for different $`N`$ values and for $`A/\mathrm{}=2000s^1`$, $`\sigma =5\mu m`$. The same distributions, but for $`\omega _{}=2\pi \times 600s^1`$, are given in Fig. 6(b).

Experimentally, evidence of SC against a trivial statistical mixture is obtainable by setting the system at the above parameter values and observing the coherent oscillation between the two energy minima by a measurement not destroying the coherence (e.g., a phase contrast technique). The condensate can be prepared in the ground state, because the system condensates naturally in this state. Then one measures the condensate barycenter by the above noninvasive technique at the initial time. This measurement collapses the distribution $`P(m)`$ onto one peak. If the energy transfer of the measurement is not too large, we expect that only the first excited state is populated and hence we have the superposition of $`\varphi _0(m)`$ and $`\varphi _1(m)`$ \[see Fig. 7(a)\]. At the following times the distribution oscillates between the two peaks. The symmetry of Hamiltonian (10) with respect to the transformation $`\widehat{a}\widehat{a}`$ shows that $`\varphi _0(m)`$ and $`\varphi _1(m)`$ are symmetric or antisymmetric; by a suitable choice of the phase, both $`\varphi _0`$ and $`\varphi _1`$ can be real functions. We consider a state preparation such that at $`t=0`$ the state $`\varphi _0(m)+\varphi _1(m)`$ has only one probability peak \[see Fig. 7(a), dashed line\]. After a quarter period corresponding to the frequency separation $`\omega _f`$ between ground and first excited states, the superposition will be $`\varphi _0(m)i\varphi _1(m)`$ and the corresponding probability is the two-peak aspect (solid line),

$$Q(m,t=\pi /\omega _f)=\varphi _0(m)^2+\varphi _1(m)^2.$$
(18)
At time $`t=\pi /\omega _f`$ the only peak is that missing at time $`t=0`$, thus there is a coherent oscillation between the two states. Detecting such an oscillation would provide evidence of an SC at an intermediate time when both peaks are present.

The intermediate distribution, given by Eq. (18), may display two peaks even when the ground state has only one \[see Fig. 7(b)\]. This occurs because the distribution $`\varphi _1(m)^2`$ has always two peaks.

Notice that if we had $`N`$ loosely coupled or independent atoms ($`Ng`$ sufficiently low or even zero) the superposition of ground and first excited state would have a single peak, oscillating as a coherent state inside a harmonic potential. This would by no means be a SC. On the contrary, we have shown that for $`Ng`$ sufficiently high, we have a two peak distribution with the two partial barycenters at nearly constant positions. During the evolution, the two peak amplitudes oscillate, that is, the probability to find the system in either state oscillates.

In conclusion, we have found that the oscillation frequency between the two states of the SC is $`\omega _f50s^1`$ \[see Fig. 2(a) for $`NN_i5`$\]. In order to neglect thermal activation, the system should be cooled at a temperature of around $`1nK`$. To cool at $`1nK`$ is within the reach of present technologies, even though such a low temperature has not yet been reported. On the other hand, the collective degrees of freedom for which there is quantum coherence may be weakly coupled with the other modes of the condensate, which act as a thermal bath. Hence, in such a case it might not be necessary to cool at $`1nK`$ the whole condensate, but just the involved degrees of freedom.

## IV Losses and observability of MQC

Besides the finite temperatures effects, we must also account for atom losses from the condensate.

The entanglement between condensate and lost atoms will imply a decoherence of the superposition. Near the threshold for symmetry breaking, the two macroscopic wave functions corresponding to the energy minima of the classical system are almost coincident (see Fig. 1, for $`N=500`$), thus the loss of one or a few atoms does not allow us yet to discriminate between the two alternatives, and thus decoherence effects are almost negligible.

We perform a simple calculation of the decoherence effect due to atom loss for the condensate. In the Hartree-Fock approximation all the atoms are in the same state, given by the GP wave function. Hence, if $`\psi _0(\stackrel{}{x})`$ is one of the two asymmetrical states with minimal energy, the wavefunction for the whole system of $`N`$ atoms is given by

$$\psi _l(\stackrel{}{x}_1,\stackrel{}{x}_2,\mathrm{},\stackrel{}{x}_N)=\underset{i=1}{\overset{N}{}}\psi _0(\stackrel{}{x}_i).$$
(19)
But the GP yields another minimum $`\psi _0(\stackrel{}{x})`$, for which the overall wave function is given by

$$\psi _r(\stackrel{}{x}_1,\stackrel{}{x}_2,\mathrm{},\stackrel{}{x}_N)=\underset{i=1}{\overset{N}{}}\psi _0(\stackrel{}{x}_i).$$
(20)
The coherent superposition of the two alternatives is given by

$$\psi _s=K(\psi _l+\psi _r)$$
(21)
where $`K`$ is a normalization factor. If $`\psi _{l,r}`$ are orthogonal then $`K=\frac{1}{\sqrt{2}}`$.

The density operator is

$$\widehat{\rho }=\frac{1}{2}(|\psi _l+|\psi _r)(\psi _l|+\psi _r|)$$
(22)
The coherence between the two alternatives is given by

$$C=2Tr[|\psi _l\psi _r|\widehat{\rho }]=1.$$
(23)
If we assume that the atoms escaping the condensate transfer a negligible energy to the trapped atoms, then the residual coherence after the loss of $`M`$ atoms is

$$\stackrel{~}{C}=2Tr[|\stackrel{~}{\psi }_l\stackrel{~}{\psi }_r|\widehat{\rho }],$$
(24)
where the vectors $`|\stackrel{~}{\psi }_{l,r}>`$ refer to the $`NM`$ atoms still in the condensate. The functions $`\psi _{l,r}`$ can be written as

$$\psi _{l,r}=\stackrel{~}{\psi }_{l,r}\underset{i=1}{\overset{M}{}}\psi _0(\pm \stackrel{}{x}_i),$$
(25)
where $`i=1,2,\mathrm{},M`$ are the lost atoms. We rewrite $`\widehat{\rho }`$ by using this expansion, and find, provided that $`\stackrel{~}{\psi }_l|\stackrel{~}{\psi }_r=0`$,

$$\stackrel{~}{C}=\left(d^3x\psi _0(\stackrel{}{x})\psi _0(\stackrel{}{x})\right)^M(I_o)^M.$$
(26)
Let us call $`ϵ=1I_o`$. If $`ϵ1`$ then the coherence is given by

$$\stackrel{~}{C}e^{ϵM}.$$
(27)
The quantity $`N_d=1/ϵ`$ provides the number of atoms that must be lost in order to reduce the coherence by $`1/e`$. Figure 8 shows how $`N_d`$ scales with the total number of atoms $`N`$ for $`A=2000s^1`$ and $`\sigma =5\mu m`$.

Notice that near the threshold value ($`N=643`$) $`N_d`$ is rather large, and hence, coherence is more robust with respect to atom losses. Far threshold $`N_d`$ has a weak dependence on $`N`$.

We now evaluate the loss rate of the condensate. The relevant processes are two-body inelastic and three-body collisional decays. The total loss rate of the two decays are Dodd

$$R(N)=\alpha N^2d^3r|\psi _0(r)|^4+LN^3d^3r|\psi _0(\stackrel{}{x})|^6$$
(28)
where $`\alpha `$ is the two-body dipolar loss rate coefficient and $`L`$ is the three-body recombination loss rate coefficient. We use the following values: $`\alpha =1.2\times 10^{14}cm^3s^1`$ Dodd , $`L=2.6\times 10^{28}cm^6s^1`$ Moe .

We report in Fig. 9 $`R(N)`$ for $`A=2000s^1`$ and $`\sigma =5\mu m`$. From Figs. 2(b),8, and 9 we can infer that near threshold atom losses have a small effect on coherence.

For instance, for $`N=650`$ the tunneling rate is about $`25s^1`$ \[see Figs. 2(b) and 3(b)\], thus the number of atoms lost during an oscillation period ($`0.25s^1`$) is around $`3.5`$ (see Fig. 9); but Fig. 8 shows that such a loss is negligible for coherence. For $`N=655`$ the tunneling rate is about $`10s^1`$, thus the number of atoms lost during an oscillation period is about $`9`$, still below $`N_d30`$.

Besides decoherence, the atom loss yields a shift in the tunneling frequencies. As shown in Figs. 2-3, the frequencies vary a lot as we reduce the number of condensed atoms. For $`N=650`$ the lost atoms are around $`3.5`$ during an oscillation period, thus the oscillation frequency changes by $`30\%`$.

In order to reduce the effect of losses we lower the threshold value by increasing the frequency $`\omega _{}`$. In Figs. 10 and 11 we report on $`N_d`$ and $`R(N)`$ for $`\omega _{}=2\pi \times 600s^1`$; the longitudinal frequency and the standing wave have been kept unchanged. The tunneling frequencies and probabilities $`P(m)`$ are given in Figs. 4(b) and 6(b), respectively. We notice that here the loss rate is much smaller, and hence the frequency shift is smaller, whereas the tunneling frequencies are slightly larger. A further reduction of the threshold value should improve the situation.

## V Conclusion

We have solved numerically the Gross-Pitaevskii equation for a $`{}_{}{}^{7}Li`$ condensate in a double-well potential and we have shown that the spatial density undergoes a symmetry breaking for a suitable number of condensed atoms. The classical asymmetrical wave function is used to select two modes. With their quantization we have evaluated the tunneling rate for some parameters of the system and we have found that it is within the reach of present laboratory technologies. Thermal effects could be lower than expected because of a decoupling between the collective variables and the thermal bath. To prove the presence of two alternatives we have introduced an appropriate observable. We have then evaluated the effects of loss and we have showed that the decoherence is negligible. In a forthcoming work we discuss the conditions for bistability in the case of mutual repulsion (positive scattering length).


format: text

text:
# Nonlinear Particle Acceleration in Relativistic Shocks

## 1 Introduction

Most collisionless shocks in astrophysics are nonrelativistic, i.e., the flow speed of the unshocked plasma in the reference frame at rest with the shock, $`u_0`$, is much less than the speed of light, $`c`$. Particle acceleration in such shocks has been studied extensively in both the linear and nonlinear regimes, and we refer the reader to reviews which describe the basic features of the shocks, the energetic particles they produce, and the many applications (e.g., Drury, 1983; Blandford & Eichler, 1987; Jones & Ellison, 1991). Relativistic shocks, where the flow speed Lorentz factor $`\gamma _0=[1(u_0/c)^2]^{1/2}`$ is greater than a few, are likely to be much less common than nonrelativistic ones, but may occur in extreme objects such as pulsar winds, hot spots in radio galaxies, and gamma-ray bursts (GRBs). Largely motivated by the application to GRBs, relativistic shocks have recently received considerable attention by a number of researchers (e.g., Bednarz & Ostrowski, 1996; Kirk et al., 2000; Achterberg et al., 2001). However, except for some preliminary work done over a decade ago (Schneider & Kirk, 1987; Ellison, 1991a, b), current descriptions of relativistic shocks undergoing first-order Fermi acceleration are test particle approximations that do not include the backreaction of the accelerated particles on the shock structure. This may be a serious limitation of relativistic shock theory in applications, such as GRBs, where high particle acceleration efficiencies are often assumed. Here, we present results for nonlinear acceleration where the backreaction of the accelerated population on the relativistic shock structure is included self-consistently.

In collisionless shocks, charged particles gain energy by scattering back and forth between the converging upstream and downstream plasmas. This basic physical process, called diffusive or first-order Fermi shock acceleration, is the same in relativistic and nonrelativistic shocks. Differences in the mathematical description and outcome of the process occur, however, because energetic particle distributions are nearly isotropic in the shock reference frame in nonrelativistic shocks (where $`vu_0`$; $`v`$ is the individual particle speed), but are highly anisotropic in relativistic shocks (since $`vu_0c`$) (e.g., Peacock, 1981). The description of particle diffusion and energy gain is far more difficult when $`\gamma _01`$ because the diffusion approximation, which requires nearly isotropic distribution functions, cannot be made. Because of this, Monte Carlo simulations, where particle scattering and transport are treated explicitly, and which, in effect, solve the Boltzmann equation with collective scattering (e.g., Ellison & Eichler, 1984; Kirk & Schneider, 1987b; Ellison, Jones, & Reynolds, 1990; Ellison & Reynolds, 1991; Ostrowski, 1991; Bednarz & Ostrowski, 1996), offer advantages over analytic methods. This is true in the test-particle approximation, where analytic results exist, but is even more important for nonlinear relativistic shocks.

In nonrelativistic shocks, for $`vu_0`$, a diffusion-convection equation can be solved directly for infinite, plane shocks (e.g., Axford, Leer, & Skadron, 1977; Blandford & Ostriker, 1978), yielding the well-known result

$$f(p)d^3pp^\sigma d^3p\text{with}\sigma =3r/(r1),$$
(1)
where $`r`$ is the shock compression ratio, $`p`$ is the momentum, and $`f(p)d^3p`$ is the number density of particles in $`d^3p`$. Eq. (1) is a steady-state, test-particle result with an undetermined normalization, but the spectral index, $`\sigma `$, in this limit is independent of the shock speed, $`u_0`$, or any details of the scattering process as long as there is enough scattering to maintain isotropy in the local frame. To obtain an absolute injection efficiency, or to self-consistently describe the nonlinear backreaction of accelerated particles on the shock structure (at least when the seed particles for acceleration are not fully relativistic to begin with), techniques which do not require $`vu_0`$ must be used. Furthermore, for particles that do not obey $`vu_0`$ additional assumptions must be made for how these particles interact with the background magnetic waves and/or turbulence, i.e., the so-called “injection problem” must be considered (see, for example, Jones & Ellison, 1991; Malkov, 1998). The Monte Carlo techniques we describe here make the simple assumption that all particles, regardless of energy, interact in the same way, i.e., all particles scatter elastically and isotropically in the local plasma frame with a mean free path proportional to their gyroradius. These techniques and assumptions have been used to calculate nonlinear effects in nonrelativistic collisionless shocks for a number of years with good success comparing model results to spacecraft observations (e.g., Ellison & Eichler, 1984; Ellison, Möbius, & Paschmann, 1990; Ellison, Jones, & Baring, 1999).

Early work on relativistic shocks was mostly analytical in the test particle approximation (e.g., Blandford & McKee, 1976; Peacock, 1981; Kirk & Schneider, 1987a; Heavens & Drury, 1988), although the analytical work of Schneider & Kirk (1987) explored modified shocks. Test-particle Monte Carlo techniques for relativistic shocks were developed by Kirk & Schneider (1987b) and Ellison, Jones, & Reynolds (1990) for parallel, steady-state shocks, i.e., those where the shock normal is parallel to the upstream magnetic field, and extended to include oblique magnetic fields by Ostrowski (1991). Some preliminary work on modified relativistic shocks using Monte Carlo techniques was done by Ellison (1991a, b).

The most important results from the theory of test-particle acceleration in ultrarelativistic shocks are: (i) regardless of the state of the unshocked plasma, some particles can pick up large amounts of energy $`\mathrm{\Delta }E\gamma _0^2mc^2`$ in their first shock crossing cycle (Vietri, 1995), but will receive much smaller energy boosts ($`<E_f/E_i>2`$) for subsequent crossing cycles (e.g., Gallant & Achterberg, 1999; Achterberg et al., 2001)111$`E_i`$ ($`E_f`$) is the particle energy at the start (end) of an upstream to downstream to upstream (or a downstream to upstream to downstream) shock crossing cycle.; (ii) the shock compression ratio, defined as $`ru_0/u_2`$, tends to $`3`$ as $`u_0c`$ (e.g., Peacock, 1981; Kirk, 1988), where $`u_2`$ is the flow speed of the shocked plasma measured in the shock frame;222Note that the density ratio across the relativistic shock $`\rho _2/\rho _0u_0/u_2`$, in contrast with nonrelativistic shocks, because the Lorentz factors associated with the relativistic flows modify the particle flux jump condition. Here and elsewhere we use the subscript 0 (2) to indicate far upstream (downstream) values. and (iii) a so-called ‘universal’ spectral index, $`\sigma 4.24.3`$ (in Eq. 1) exists in the limits of $`\gamma _01`$ and $`\delta \theta 1`$, where $`\delta \theta `$ is the change in direction a particle’s momentum vector makes at each pitch angle scattering (e.g., Bednarz & Ostrowski, 1998; Achterberg et al., 2001).

We find that these results are modified in mildly relativistic shocks, even in the test-particle approximation, and in fully relativistic shocks (at least for $`\gamma _010`$) when the backreaction of the accelerated particles is treated self-consistently, which causes the shock to smooth and the compression ratio to change from test-particle values. In mildly relativistic shocks, $`f(p)`$ remains a power law in the test-particle approximation but both $`r`$ and $`\sigma `$ depend on the shock Lorentz factor, $`\gamma _0`$. When efficient particle acceleration occurs in mildly relativistic shocks (i.e., $`\gamma _03`$), large increases in $`r`$ can result and a power law is no longer a good approximation to the spectral shape. In these cases, we determine the compression ratio by balancing the momentum and energy fluxes across the shock with the Monte Carlo simulation. For larger Lorentz factors, accelerated particles smooth the shock structure just as they do in slower shocks, but $`r`$ approaches 3 as $`\gamma _0`$ increases. In general, efficient particle acceleration results in spectra very different from the so-called ‘universal’ power law found in the test-particle approximation unless $`\gamma _010`$.

## 2 Monte Carlo Model

The techniques we use are essentially identical to those described in Ellison, Baring, & Jones (1996) and Ellison, Jones, & Baring (1999). The differences are that the code has been made fully relativistic and only results for parallel shocks with pitch-angle diffusion are presented here. The code is steady-state, includes a uniform magnetic field, and moves particles in helical orbits. We assume the Alfvén Mach number is large, i.e., we neglect any effects from Alfvén wave heating in the upstream precursor. This also means we neglect the second-order acceleration of particles scattering between oppositely propagating Alfvén waves. Such an effect in relativistic plasmas with strong magnetic fields is proposed for nonlinear particle acceleration in GRBs by Pelletier (1999) (see also Pelletier & Marcowith, 1998).

The pitch angle diffusion is performed as described in (Ellison, Jones, & Reynolds, 1990). That is, after a small increment of time, $`\delta t`$, a particles’ momentum vector, $`𝐩`$, undergoes a small change in direction, $`\delta \theta `$. If the particle originally had a pitch angle, $`\theta `$ (measured relative to the shock normal direction), it will have a new pitch angle $`\theta ^{}`$ such that

$$\mathrm{cos}\theta ^{}=\mathrm{cos}\theta \mathrm{cos}\delta \theta +\sqrt{1\mathrm{cos}^2\theta }\mathrm{sin}\delta \theta \mathrm{cos}\varphi ,$$
(2)
where $`\varphi `$ is the azimuth angle measured with respect to the original momentum direction. All angles are measured in the local plasma frame. If $`\delta \theta `$ is chosen randomly from a uniform distribution between 0 and $`\delta \theta _{\mathrm{max}}`$ and $`\varphi `$ is chosen from a uniform distribution between 0 and $`2\pi `$, the tip of the momentum vector will perform a random walk on a sphere of radius $`p`$. As shown by Ellison, Jones, & Reynolds (1990), the angle $`\delta \theta _{\mathrm{max}}`$ is determined by

$$\delta \theta _{\mathrm{max}}=\left(6\delta t/t_c\right)^{1/2}=\left(12\pi /N\right)^{1/2},$$
(3)
where $`N=\tau _g/\delta t1`$ is the number of gyro-segments, $`\delta t`$, dividing a gyro-period $`\tau _g=2\pi r_g/v`$. The time $`t_c`$ is a “turn around” time defined as $`t_c=\lambda /v`$, where $`\lambda `$ is the particle mean free path. The mean free path is taken to be proportional to the gyroradius $`r_g=pc/(QeB)`$ ($`e`$ is the electronic charge, $`Q`$ is the ionic charge number, and $`B`$ is the local uniform magnetic field), i.e., $`\lambda =\eta r_g`$, where $`\eta `$ determines the strength of scattering. In all of the examples given here we set $`\eta =1`$, or, in other words, the strong scattering Bohm limit is assumed.

For a downstream particle to return upstream, it’s velocity vector must be directed within a cone with opening angle $`\theta _2`$ such that $`|v_2\mathrm{cos}\theta _2|>u_2`$, where $`v_2`$ and $`\theta _2`$ are measured in the downstream frame and $`\theta _2=0^{}`$ is in the $`x`$-direction, i.e., along the shock normal direction. For fully relativistic shocks with $`v_2c`$ and $`u_2c/3`$, $`\mathrm{cos}\theta _21/3`$ for a downstream particle to cross the shock into the upstream region. When the particle enters the upstream region it must satisfy essentially the same constraint, i.e., $`|v_0\mathrm{cos}\theta _0|>u_0`$, where now $`v_0`$ and $`\theta _0`$ are measured in the upstream frame.333In the test-particle approximation, $`u_0`$ is just the shock speed. In nonlinear shocks, the flow speed just upstream from the subshock at $`x=0`$ will be less than the far upstream shock speed, $`u_0`$, as measured in the shock reference frame. Since both the particle and shock have high Lorentz factors, we can write

$$\mathrm{cos}\theta _0=u_0/v_0=\frac{\left(11/\gamma _0^2\right)^{1/2}}{\left(11/\gamma _v^2\right)^{1/2}}\frac{11/(2\gamma _0^2)}{11/(2\gamma _v^2)},$$
(4)
where $`\gamma _v[1(v_0/c)^2]^{1/2}`$ is the particle Lorentz factor. Since $`\mathrm{cos}\theta _01\theta _0^2/2`$ for small $`\theta _0`$, we have

$$\theta _0^2\frac{1}{\gamma _0^2}\frac{1}{\gamma _v^2}.$$
(5)
For ultrarelativistic particles with $`\gamma _v\gamma _0`$, $`\theta _01/\gamma _0`$ (e.g., Gallant & Achterberg, 1999), but $`\theta _0`$ can be much smaller for mildly relativistic particles.

In order to re-cross into the downstream region, particles must scatter out of the upstream cone defined by $`\theta _0`$ and Achterberg et al. (2001) show that most particles are only able to change the angle they make with the upstream directed shock normal by $`|\delta \theta |\theta _0`$ before being sweep back downstream, making the distribution of shock crossing particles highly anisotropic. Therefore, if the shock Lorentz factor $`\gamma _01`$, a larger fraction of particles re-cross the shock into the downstream region with highly oblique angles (as measured in the shock frame) compared to lower speed shocks (see Fig. 2 discussed below). Particles crossing at such oblique angles receive smaller energy gains than would be the case for an isotropic pitch angle distribution and Achterberg et al. (2001) go on to show that $`<E_f/E_i>2`$ for a shock crossing cycle (after the first one).

Considering these constraints, we require $`\delta \theta _{\mathrm{max}}<\theta _0`$, or

$$N>N_{\mathrm{min}}=12\pi \gamma _0^2,$$
(6)
and find (as shown in Fig. 1) that the power law spectral index, $`\sigma `$, asymptotically approaches a maximum value as $`N`$ is increased. If $`N`$ is less than the value required for convergence (and the gyro-segments are too large), the distribution will be flatter than produced with the convergent value of $`N`$ because more particles are able to cross from upstream to downstream with $`\theta _{\mathrm{sk}}90^{}`$ and receive unrealistically large energy boosts. This effect has long been known from the comparison of pitch-angle diffusion to large-angle scattering in relativistic shocks (e.g., Kirk & Schneider, 1987b; Ellison, Jones, & Reynolds, 1990). For all of the examples reported on here, $`N`$ is chosen large enough so it makes no difference if $`\delta \theta `$ is chosen uniformly between $`0^{}`$ and $`\delta \theta _{\mathrm{max}}`$ or if $`\mathrm{cos}\delta \theta `$ is chosen uniformly between $`\mathrm{cos}\delta \theta _{\mathrm{max}}`$ and 1. Fig. 1 shows how our results depend on $`N`$ for shock speeds ranging from fully nonrelativistic to fully relativistic. In all cases, as $`N`$ is increased the spectral index approaches a maximum and for $`\gamma _07`$ we obtain the well known result $`\sigma =4.2`$$`4.3`$. The fact that the computation time for the Monte Carlo simulation scales as $`N`$ and $`N\gamma _0^2`$ places limits on modeling ultrarelativistic shocks with this technique.

In Fig. 2 we compare pitch-angle distributions (measured in the shock reference frame) of particles crossing the shock. The curves are normalized such that the area under each curve equals one and the $`x`$-component of particle momentum in the shock frame, $`p_x`$, is positive when directed downstream to the right (see Fig. 8 for the shock geometry). Here, $`p_t`$ is the magnitude of the total particle momentum also measured in the shock frame. In the top panel, we compare an unmodified (UM) and nonlinear (NL) mildly relativistic shock ($`\gamma _0=2`$) with a nonrelativistic one ($`u_0=5000`$ km s-1). Particles crossing the $`\gamma _0=2`$ shock are highly anisotropic with $`p_x/p_t`$ strongly peaked near $`0.35`$. In the nonrelativistic shock, the particles are nearly isotropic except for a slight flux-weighting effect. There is little difference in the distributions between the UM and NL shocks. In the bottom panel, we show the pitch-angle distributions for UM and NL shocks with $`\gamma _0=10`$. While the distributions are somewhat more sharply peaked, they are quite similar to those for $`\gamma _0=2`$ and show relatively small variations between the UM and NL shocks.

The main difference between our present code and an earlier code used by Ellison, Jones, & Reynolds (1990) to model test-particle relativistic shocks is that the previous code used a guiding center approximation with an emphasis on large-angle scattering rather than the more explicit orbit calculation of pitch-angle diffusion used here. Other than the far greater range in $`\gamma _0`$ and the nonlinear results we now present, the work of Ellison, Jones, & Reynolds (1990) is consistent with the work presented here.

The particle transport is performed as follows. Particles of some momentum, $`p_{\mathrm{pf}}`$, (measured in the local plasma frame) are injected far upstream from the shock and pitch-angle diffuse and convect until they cross a grid zone boundary, i.e., a dividing plane in the simulation between regions with different bulk flow speeds (see Ellison, Baring, & Jones, 1996, for a full discussion). For unmodified shocks there is only one grid zone boundary which divides the upstream and downstream regions, but for nonlinear shocks the bulk flow speed changes in small steps, each separated by a boundary, from $`u_0`$ far upstream to $`u_2`$ downstream. In the unmodified case, the shock thickness is essentially zero (i.e., shorter than the distance a particle diffuses in $`\delta t`$) but in the nonlinear, modified case, the shock precursor extends over the entire region of varying bulk flow speeds and a small scale “subshock” (at position $`x=0`$ in our simulation) exists where most of the entropy production occurs. When a particle crosses a grid zone boundary, $`p_{\mathrm{pf}}`$ is transformed to the new local frame moving with a new speed relative to the subshock and the particle continues to scatter and convect. Each particle is followed until it leaves the system in one of three ways. It can convect far downstream and not return to the subshock, it can obtain a momentum greater than some $`p_{\mathrm{max}}`$ and be removed, or it can diffuse far enough upstream to cross an upstream free escape boundary (FEB) and be removed. Both $`p_{\mathrm{max}}`$ and the position of the FEB are free parameters in our model (see Berezhko & Völk, 1997, for a discussion of the self-consistent determination of the maximum particle energy in nonrelativistic supernova remnant shocks). For our nonlinear calculations, we iterate the shock structure and compression ratio until the number, momentum, and energy fluxes are conserved across the shock. This procedure has been detailed many times for nonrelativistic shocks (see Ellison, Baring, & Jones, 1996, and references therein) and the modifications required for relativistic shocks were given in Ellison & Reynolds (1991).

To avoid excessive computation, we use a probability of return calculation as described in detail in Ellison, Jones, & Reynolds (1990). That is, we use the standard expression obtained by Peacock (1981),

$$P_R=\left(\frac{vu_2}{v+u_2}\right)^2,$$
(7)
to determine the probability, $`P_R`$, that a particle, having crossed a particular point in the uniform downstream flow, will return back across that point. Eq. (7) is fully relativistic and independent of the diffusive properties of the particles as long as they are isotropic in the $`u_2`$ frame. We ensure this isotropy by only applying Eq. (7) once a particle has diffused several mean free paths downstream from the shock. For ultrarelativistic shocks with $`vu_0c`$ and $`r3`$, Eq. (7) gives $`P_R0.25`$.

For clarity, we note that this is not the same probability, $`𝒫_{\mathrm{ret}}`$, that is used by Achterberg et al. (2001) to determine the test-particle power law index, i.e.,

$$s=1+\frac{\mathrm{ln}(1/𝒫_{\mathrm{ret}})}{\mathrm{ln}<E_f/E_i>},$$
(8)
where $`N(E)E^s`$ and, for fully relativistic particles, $`s=\sigma 2`$. The quantity $`<E_f/E_i>`$ is the average energy ratio for a particle undergoing a shock crossing cycle. Although not explicitly stated in Achterberg et al. (2001), it is clear from the context that $`𝒫_{\mathrm{ret}}`$ is calculated just behind the shock where the particle distribution is highly anisotropic. In this case, particles that have just crossed from upstream to downstream will be more likely to recross back into the upstream region than indicated by Eq. (7) because their pitch angles are more likely to be highly oblique relative to the shock normal than in the isotropic distributions further downstream (compare the solid or dashed curves to the dotted curve in the top panel of Fig. 2 discussed below). For $`\gamma _0=10`$, Achterberg et al. (2001) find $`𝒫_{\mathrm{ret}}=0.435\pm 0.005`$ and $`<E_f/E_i>=1.97\pm 0.01`$ giving the standard result $`s=2.230\pm 0.012`$. As shown in Fig. 1, our unmodified results are consistent with this spectral index for $`\gamma _07`$.

## 3 Test-Particle Results

In Fig. 3 we show particle distributions for unmodified shocks with speeds ranging from fully nonrelativistic ($`u_0=5000`$ km s-1) to mildly relativistic ($`\beta _0=u_0/c=0.5`$) to fully relativistic ($`\gamma _0=10`$). The nonrelativistic distribution matches the standard test-particle Fermi result of $`\sigma =4`$ for $`r=4`$ and the fully relativistic result is consistent with the well-known limit of $`\sigma \mathrm{\hspace{0.17em}4.2}4.3`$ as $`\gamma _0\mathrm{}`$ for $`r=3`$. In the trans-relativistic regime, both the compression ratio and the spectral index vary with $`u_0`$. For the $`u_0=0.5c`$ distribution shown in Fig. 3, we have determined the compression ratio by balancing the mass, momentum, and energy fluxes across the shock in the test-particle limit, i.e., by ignoring any effects from the accelerated particles. This technique is described in detail in Ellison & Reynolds (1991). We find $`r=3.8\pm 0.1`$ and the spectral index is $`\sigma 4.01`$, slightly flatter than $`3r/(r1)4.07`$, the nonrelativistic result.

Fig. 4 shows the compression ratio as a function of $`\beta _0\gamma _0`$ (solid curve), still ignoring the effects of accelerated particles. The compression ratio is determined self-consistently by balancing the momentum and energy fluxes across the shock with no restriction on the shocked or unshocked adiabatic index. As expected, $`r`$ decreases smoothly from $`4`$ for fully nonrelativistic shocks to $`3`$ for fully relativistic shocks. The power law index, $`\sigma `$, is also shown (dashed curve) and this varies slowly from $`\sigma =4`$ to $`\sigma 4.23`$ between the two extremes. For comparison, we show $`r`$ (solid dots) for shocks undergoing efficient particle acceleration. For these points, all shock parameters except $`\beta _0\gamma _0`$ and $`p_{\mathrm{max}}`$ are kept constant. For the nonlinear shocks, the maximum cutoff momentum is set to $`p_{\mathrm{max}}=10^4\beta _0\gamma _0m_pc`$ in all cases.444The value of $`p_{\mathrm{max}}`$ can have a large influence on the shock characteristics at low $`\beta _0\gamma _0`$ because $`r`$ is large enough (and $`\sigma `$ small enough) that particles escaping at $`p_{\mathrm{max}}`$ are dynamically important. When $`\beta _0\gamma _0`$ becomes large enough so that $`r`$ drops below $`4`$, $`p_{\mathrm{max}}`$ is no longer an important parameter for the shock structure. We discuss the nonlinear results in detail below, but here only emphasize that $`r`$ in nonlinear shocks will be larger than the test-particle value for $`\gamma _010`$. The test-particle results shown in Fig. 4 are in close agreement with recent analytic results of Kirk et al. (2000) and Gallant (2002). To our knowledge, no analytic results exist for nonlinear, relativistic shocks.

In Fig. 5 we show the average ratios of momenta (measured in the local plasma frame) for particles executing upstream to downstream to upstream cycles across the shock, $`<p_f/p_i>_{\mathrm{u}\mathrm{d}\mathrm{u}}`$, and downstream to upstream to downstream cycles, $`<p_f/p_i>_{\mathrm{d}\mathrm{u}\mathrm{d}}`$. These results are for $`\gamma _0=10`$ ($`r=3`$) and show a slight momentum dependence at low momenta but converge to $`<p_f/p_i>=2.2\pm 0.1`$ at high momenta. This value is close to $`<E_f/E_i>=1.97\pm 0.01`$ reported by Achterberg et al. (2001) for $`\gamma _0=10`$. As mentioned above, in the first upstream to downstream to upstream cycle, particles achieve a large boost in momentum as indicated in the figure.

The difference between our value of $`<p_f/p_i>`$ at large $`p_i`$ and that of Achterberg et al. (2001) is greater than the uncertainties and probably stems from the different assumptions made in the simulations. As discussed above, $`<p_f/p_i>`$ depends critically on the average angle a particle makes when crossing the shock. While the large majority of particles in our simulations gain energy by crossing from downstream to upstream and then immediately (within a few $`\delta t`$’s) re-crossing back into the downstream region at oblique angles, a few manage to diffuse farther upstream (see Fig. 7 below). When these particles re-cross the shock into the downstream region, they can do so at flatter angles and receive larger energy gains. We speculate that differences in how these few particles are treated in the simulations might produce the differences in $`<p_f/p_i>`$.

As an illustration of how particles interact with nonrelativistic and relativistic unmodified shocks, we shown trajectories for two individual particles in Figs. 6 and 7. The lower panel in each figure shows a trace of the particle trajectory and the upper panels show the particle momentum, always measured in the local plasma frame. For the nonrelativistic shock (Fig. 6), the speed of the particle is far greater than the shock speed and it diffuses easily on both sides of the shock. When it crosses $`x=0`$, it does so nearly isotropically (except flux weighting makes crossings with flat trajectories slightly more likely) and essentially always gains momentum. The momentum gain in a single shock crossing is small, but a particle can stay in the system for many crossings.

When the shock speed, $`u_0`$, is close to $`c`$, the particle will be convected downstream much more rapidly than in nonrelativistic shocks and few particles will be able to cross the shock many times. However, downstream particles that do manage to cross the shock into the upstream region do so with much flatter trajectories, as discussed above, and can receive large momentum boosts in a single shock crossing due to the shock’s Lorentz factor (note the logarithmic scale in the top panel of Fig. 7). In a typical shock crossing cycle, downstream particles gain momentum when they cross into the upstream region (see positions labeled $`a`$ and $`b`$ in Fig. 7), lose momentum when they cross back downstream because they cross with oblique pitch angles, but end up with a net momentum boost. However, as shown by the position labeled $`c`$, it is possible for a particle to diffuse farther upstream before being convected back to the shock. In this case, it can cross the shock with a flat trajectory and gain momentum upon entering the downstream region. If the acceleration is efficient, the few particles that diffuse far upstream carry enough pressure to produce the shock smoothing we discuss next.

## 4 Non-Linear Results

As noted by Achterberg et al. (2001), the large energy boost particles receive in their initial crossing of the shock provides a natural injection process for further acceleration and suggests that relativistic shocks may be efficient accelerators. However, just as with nonrelativistic shocks, efficient acceleration limits the use of test-particle approximations and requires that the nonlinear backreaction of the accelerated particles be treated self-consistently (e.g., Jones & Ellison, 1991). These nonlinear effects will result in a smoothing of the shock and a change in the overall shock compression ratio, just as they do in nonrelativistic shocks. The differences between test-particle and nonlinear results, for the parameter ranges we have investigated, are large enough to produce spectra noticeably different from the often quoted $`N(E)E^{2.3}`$ and to influence applications to GRB models where high particle acceleration efficiencies are assumed. We illustrate this with two examples, one mildly relativistic ($`\gamma _0=1.4`$) and one more fully relativistic ($`\gamma _0=10`$). The far upstream conditions have relatively little influence on our results as long as $`\gamma _p\gamma _0`$, where $`\gamma _p`$ is the plasma frame Lorentz factor of the far upstream, injected particles. For concreteness, in both examples we take the far upstream plasma to be a thermal distribution of protons at a temperature of $`10^6`$ K.

### 4.1 Mildly relativistic shock: $`\gamma _0=1.4`$

Fig. 8 shows unmodified and nonlinear shock structures for $`\gamma _0=1.4`$. The top panel shows $`\gamma _u(x)u(x)`$, where $`\gamma _u(x)=\{1[u(x)/c]^2\}^{1/2}`$, the middle panel is momentum flux, and the bottom panel is the energy flux, all scaled to far upstream values. All curves are plotted versus $`x`$, where $`x=0`$ is the position of the sharp subshock. A logarithmic scale is used for $`x<10r_{g,0}`$ and a linear scale is used for $`x>10r_{g,0}`$, where $`r_{g,0}m_pu_0/(eB)`$.555In these parallel shock simulations where we ignore Alfvén wave production, the magnetic field, $`B`$, has no other effect than setting arbitrary length and time scales. In each panel, the solid curve is from an unmodified shock with $`r3.6`$ (see Fig. 4), while the dashed curve is the momentum and energy flux conserving result.

For our pitch angle diffusion model, where particles interact elastically and isotropically in the local frame according to Eq. (2) and the discussion following it, particles are accelerated efficiently enough at the unmodified shock that the momentum and energy fluxes are not conserved and rise well above the allowed far upstream values. In order to conserve these fluxes, the shock structure must be smoothed and the overall compression ratio increased above the test-particle value. Our computational scheme calculates this compression ratio and flux conserving profile and the result is the dashed curve in the top panel with the corresponding momentum and energy fluxes in the middle and bottom panels.

The source of the non-conservation of momentum and energy is the efficient acceleration of particles by the sharp flow speed discontinuity. While the actual injection and acceleration efficiency depends on our particular pitch angle diffusion model, once our scattering assumptions are made, the kinematics determine the injection and acceleration of the particles without additional parameters.666Other input parameters, such as the Mach number and $`p_{\mathrm{max}}`$ or the position of the FEB, influence the acceleration efficiency, but these are “environmental” parameters rather than parameters needed to describe the plasma interactions. Of course it would have been possible to make assumptions which resulted in an acceleration efficiency low enough that momentum and energy are approximately conserved without significantly smoothing the shock structure or changing the compression ratio from the test-particle value. For example, we could have only allowed shocked particles above some Lorentz factor $`\gamma _{\mathrm{inj}}`$ to recross into the upstream region or arbitrarily restricted the number of particles that recrossed into the upstream to a small fraction, $`f_{\mathrm{inj}}`$, of all downstream particles. By making $`f_{\mathrm{inj}}`$ low enough or $`\gamma _{\mathrm{inj}}`$ high enough we could make the efficiency as low as we wanted. However, there are at least three reasons for not making such an assumption. The first is that restricting the acceleration efficiency requires additional parameters (i.e., $`\gamma _{\mathrm{inj}}`$ and/or $`f_{\mathrm{inj}}`$) to those needed to describe pitch angle diffusion.777Alternatively, a far more complex model of the plasma interactions can be postulated than done here, inevitably requiring additional parameters (e.g., Malkov, 1998). The second is that models with inefficient acceleration will not help explain GRBs (or other objects) that require high efficiencies. If relativistic shocks are inefficient accelerators they are not very interesting. If they are efficient, they will have a qualitative resemblance to the results we show even if, as is likely, the actual plasma processes are far more complex than the simple model we use. A third, perhaps less compelling, reason is that identical scattering assumptions as used here have been used for some time in nonrelativistic shocks and shown to match both spacecraft observations (e.g., Ellison, Möbius, & Paschmann, 1990) and hybrid plasma simulations (e.g., Ellison et al., 1993) of collisionless shocks.

The increase in compression ratio from $`r3.6`$ to $`4.5`$ shown in Fig. 8 comes about, in part, because the particles escaping at $`p_{\mathrm{max}}`$ carry away momentum and energy fluxes which make the shocked plasma more compressible. This effect is countered to some degree by the fact that the self-consistent shock produces a downstream distribution with a smaller fraction of relativistic particles than the unmodified shock so that the downstream adiabatic index $`\mathrm{\Gamma }>4/3`$. This tends to produce a smaller compression ratio. The escaping fluxes show up as a lowering of the dashed curves below the far upstream values, as shown in the bottom two panels of Fig. 8, and amount to about 3% of the far upstream values for both momentum and energy. Including the escaping fluxes, the momentum and energy fluxes are conserved to within a few percent of the far upstream values. While the changes seen here are similar to those seen and discussed for many years in efficient, nonrelativistic shock acceleration (see Berezhko & Ellison, 1999, and references therein), it must be noted that there are no known analytic expressions relating escaping fluxes and $`\mathrm{\Gamma }`$ to $`r`$ in this trans-relativistic regime. One obvious difference is that for nonrelativistic shocks, the escaping momentum flux is generally much less than the escaping energy flux (when both are measured as fractions of incoming flux) since $`\rho _ev_e^3/(\rho _0u_0^3)\rho _ev_e^2/(\rho _0u_0^2)`$ when $`v_eu_0`$ (see Ellison, 1985). Here, $`v_e`$ is the velocity of the escaping particle. In relativistic shocks, $`v_eu_0c`$ so the escaping fluxes are about equal as shown in Fig. 8.

In Fig. 9 we plot $`p^4f(p)`$ for our $`\gamma _0=1.4`$ shock. The solid curve is from an unmodified (UM) shock and the dashed curve is the nonlinear (NL) result. The shock smoothing and increase in $`r`$ produce substantial differences in the spectra even though both shocks have exactly the same input conditions. (i) The overall normalization of the NL spectrum is less, reflecting the conservation of energy flux. (ii) The NL result has the distinctive concave curvature seen in nonrelativistic shocks stemming from the fact that higher momentum particles have a longer upstream diffusion length and get accelerated more efficiently than lower momentum particles in the smooth shock. (iii) The slope at the highest momentum in the NL spectrum reflects the overall compression ratio and is flatter than the TP spectrum because $`r`$ is greater. (iv) The “thermal” peak is shifted to lower momentum in the NL result and contains a larger fraction of mildly relativistic particles than in the UM result, i.e., $`\mathrm{\Gamma }1.41`$ for the NL shock compared to $`\mathrm{\Gamma }1.36`$ for the UM shock.

### 4.2 Fully relativistic, nonlinear shock: $`\gamma _0=10`$

Fig. 10 shows results for $`\gamma _0=10`$. In all panels, the solid curves are for the unmodified shock, while the dashed curves are for the flux conserving smoothed shock, both with identical input parameters. As before, the acceleration is truncated with a $`p_{\mathrm{max}}=10^4\beta _0\gamma _0m_pc9.3\times 10^{13}`$ eV/c. In Fig. 11, $`u(x)`$ and $`\gamma _u(x)`$ are plotted separately. Even though the length-scale of the shock smoothing is only a few $`r_{g,0}`$ and the change in $`r=3.03\pm 0.01`$ is small, they bring the momentum and energy fluxes into balance from values a factor of five too large in the unmodified shock. The $`1`$% difference between $`r=3.03\pm 0.01`$ and the canonical value of 3 is significant; a self-consistent solution does not exist outside of this range.

As discussed above, both the shock structure and the overall compression ratio may be modified when particle acceleration is efficient. For mildly relativistic shocks, such as our $`\gamma _0=1.4`$ example, a unique solution can be determined directly from the conserved momentum and energy fluxes. For larger $`\gamma _0`$’s however, the fluxes are less sensitive to changes in the structure and compression ratio and an additional constraint is needed to obtain a unique solution. This constraint is provided by the relativistic jump conditions for momentum and energy, i.e.,

$$\gamma _{u0}^2w_0\frac{u_0^2}{c^2}+P_0=\gamma _{u2}^2w_2\frac{u_2^2}{c^2}+P_2;$$
(9)
$$\gamma _{u0}^2w_0u_0=\gamma _{u2}^2w_2u_2,$$
(10)
where $`w=e+P`$ is the enthalpy density, $`e`$ is the total energy density, $`P`$ is the pressure, and we assume escaping fluxes are negligible. The energy density and pressure are related through a combination of the adiabatic equation of state and the conservation of energy, i.e.,

$$P=(\mathrm{\Gamma }1)(e\rho c^2),$$
(11)
where $`\rho c^2`$ is the rest mass energy density and $`\mathrm{\Gamma }`$ is, in special cases, the ratio of specific heats (e.g., Ellison & Reynolds, 1991). Dividing Eq. (9) by Eq. (10) yields

$$\frac{u_0}{c^2}+\frac{P_0}{\gamma _{u0}^2(e_0+P_0)u_0}=\frac{u_2}{c^2}+\frac{P_2}{\gamma _{u2}^2(e_2+P_2)u_2},$$
(12)
or, in terms of beta’s,

$$\beta _0+\left(\frac{P_0}{e_0+P_0}\right)\frac{1\beta _0^2}{\beta _0}=\beta _2+\left(\frac{P_2}{e_2+P_2}\right)\frac{1\beta _2^2}{\beta _2}.$$
(13)
The second term on the left hand side is small compared to $`\beta _0`$ for unshocked upstream particle temperatures less than $`10^9`$K and upstream particle densities less than $`100`$ cm-3. Neglecting this term, Eq. (13) becomes

$$\beta _0=\beta _2+\left(\frac{P_2}{e_2+P_2}\right)\frac{1\beta _2^2}{\beta _2},$$
(14)
or,

$$e_2\beta _2^2\beta _0(e_2+P_2)\beta _2+P_2=0,$$
(15)
which has the shock solution

$$\beta _2=\beta _0/r=\frac{1}{2e_2}\left[\beta _0(e_2+P_2)\sqrt{\beta _0^2(e_2+P_2)4e_2P_2}\right].$$
(16)
For $`\gamma _010`$, $`\beta _01`$ and $`re_2/P_2`$. Furthermore, if ultrarelativistic downstream particle speeds can be assumed so that $`e_2\rho _2c^2`$, Eq. (16) reduces to,

$$r\frac{1}{\mathrm{\Gamma }_21}.$$
(17)
For $`\mathrm{\Gamma }_2=4/3`$, $`r3`$, the standard ultrarelativistic result.

To obtain a unique shock solution for $`\gamma _03`$, we modify both the shock structure and the compression ratio, $`r_{\mathrm{MC}}`$, for the Monte Carlo simulation, calculate $`P_2`$ and $`e_2`$ from the resultant downstream particle distributions, and check that $`r_{\mathrm{MC}}r`$ as determined from Eq. (16). If $`r_{\mathrm{MC}}r`$, the shock structure and $`r_{\mathrm{MC}}`$ are varied until a consistent solution is found with $`r_{\mathrm{MC}}r`$. For smaller $`\gamma _0`$’s, Eq. (16) doesn’t apply because escaping fluxes become significant. In this case, however, the changes in the momentum and energy fluxes from changes in the shock structure and $`r_{\mathrm{MC}}`$ are large enough that a unique solution can be found easily, as in our $`\gamma _0=1.4`$ example.

Despite the fact that $`u(x)`$ is modified on a fairly small length scale, the resultant particle distribution function is changed substantially, as indicated in Fig. 12. In this figure, the solid curve is from the unmodified shock and the dashed curve is from the nonlinear shock, both having exactly the same input conditions. The unmodified spectrum in Fig. 12 is similar to that shown with a dashed curve in Fig. 3 only now the far upstream plasma is taken to be a thermal gas at a temperature of $`10^6`$ K rather than a delta function distribution of particles with speeds, $`v_{\mathrm{inj}}=(2E_{\mathrm{inj}}/m_p)^{1/2}`$, with $`E_{\mathrm{inj}}=1`$ MeV, as was assumed for the example in Fig. 3. The shock smoothing has caused the low energy portion of the distribution to steepen and the overall intensity to decrease, reflecting the fact that the NL spectrum conserves momentum and energy while the UM one doesn’t. The peaks in the two distributions at low momenta are also very different, with the NL spectrum having a larger fraction of slower particles than the UM one. These peaks result from the first shock crossing where all particles receive a large energy gain. In the UM case, a far greater fraction of the accelerated downstream particles are able to receive further energization by recrossing back into the upstream region than in the NL shock. The different speed distributions result in different $`\mathrm{\Gamma }`$’s and we find, by directly calculating $`\mathrm{\Gamma }`$ from the distributions, that $`\mathrm{\Gamma }_{\mathrm{UM}}1.34`$ while $`\mathrm{\Gamma }_{\mathrm{NL}}1.36`$.

### 4.3 Acceleration efficiency

The absolute acceleration efficiency can be determined in our self-consistent, nonlinear examples directly from the particle distributions. In Fig. 13 we plot, $`ϵ(>p)`$, i.e., the fraction of kinetic energy density above a momentum $`p`$ versus $`p`$ for our $`\gamma _0=1.4`$ and $`\gamma _0=10`$ examples. The fraction of kinetic energy in the quasi-thermal part of the distribution can be determined from the relatively sharp fall off of the distributions at low momenta. This occurs near $`m_pc`$ for $`\gamma _0=1.4`$ and near $`10m_pc`$ for $`\gamma _0=10`$. If we somewhat arbitrarily define the acceleration efficiency for these two examples to be $`ϵ_{1.4}(>m_pc)`$ and $`ϵ_{10}(>10m_pc)`$, respectively, we have $`ϵ_{1.4}(>m_pc)0.7`$ and $`ϵ_{10}(>10m_pc)0.6`$. Of course, the behavior of $`ϵ(>p)`$ depends on the particle spectrum from which it is derived and thus $`ϵ(>p)`$ is flatter for the $`\gamma _0=1.4`$ shock, with $`r4.5`$, than for the $`\gamma _0=10`$ shock, where $`r3.03\pm 0.01`$. Furthermore, $`ϵ(>p)`$ depends strongly on $`p_{\mathrm{max}}`$ in the $`\gamma _0=1.4`$ shock where particle escape plays an important role in determining $`r`$. The maximum momentum has little influence for $`\gamma _0=10`$ because of the steep spectrum. The fraction of kinetic energy density above $`10^3m_pc`$ is $`0.3`$ for $`\gamma _0=1.4`$ and $`0.1`$ for $`\gamma _0=10`$.

## 5 Summary and Conclusions

Particles gain energy in collisionless shocks by scattering nearly elastically off magnetic turbulence, back and forth between the converging plasmas upstream and downstream from the shock. While this basic shock acceleration physics is independent of the speed of the shock, the mathematical modeling of the process depends critically on whether or not the acceleration is efficient and whether or not particle speeds, $`v`$, are large compared to the shock speed, $`u_0`$. The Monte Carlo techniques discussed here, which do not require $`vu_0`$, are well suited for the study of relativistic shocks, and for any shock where nonlinear effects are important and the energetic particles originate as thermal particles in the unshocked plasma. Except for computational limits, these techniques allow calculations of efficient particle acceleration in shocks of any Lorentz factor.

As a check of our code, we have demonstrated that we obtain the well known, test-particle power laws in fully nonrelativistic and ultrarelativistic parallel shocks (Figs. 1 and 3). In trans-relativistic shocks, however, no such canonical results exist because the shock compression ratio, $`r`$, depends on the upstream conditions in a nonlinear fashion, even for test-particle shocks where all effects of accelerated particles are ignored (e.g., Ellison & Reynolds, 1991). We show how $`r`$, and the resulting power law spectral index, $`\sigma `$, vary through the trans-relativistic regime in Fig. 4, where $`r`$ has been determined by balancing the momentum and energy fluxes across the shock in an iterative process. These test-particle results with self-consistently determined compression ratios and adiabatic indexes are in close agreement with the recent analytic results of Kirk et al. (2000) and Gallant (2002).

Despite the fact that relativistic shock theory has concentrated almost exclusively on test-particle acceleration, it is likely that relativistic shocks are not test-particle but inject and accelerate particles efficiently. The reason is that regardless of the ambient far upstream conditions, all particles that are over taken by an ultrarelativistic shock will receive a large boosts in energy $`\mathrm{\Delta }E>\gamma _0mc^2`$ in their first shock crossing. Thus, virtually all of the particles in the downstream region of an unmodified shock are strongly relativistic with $`vc`$. The ability to overtake the shock from downstream and be further accelerated depends only on the particle speed (Eq. 7) and the presence of magnetic waves or turbulence with sufficient power in wavelengths on the order of the particle gyroradii to isotropize the downstream distributions. It is generally assumed that the necessary magnetic turbulence is self-generated and if enough turbulence is generated to scatter high momentum particles (with very low densities) that constitute a test-particle power law, there should be enough generated to isotropize lower momentum particles (which carry the bulk of the density). If acceleration can occur at all, we believe it is likely to occur efficiently making it necessary to calculate the shock structure and particle acceleration self-consistently. Furthermore, if relativistic shock theory is to be applied to gamma-ray bursts, where high conversion efficiencies are generally assumed, nonlinear effects must be calculated.

When energetic particles are generated in sufficient numbers, the conservation of momentum and energy requires that their backpressure modify the shock structure. Two basic effects occur: a precursor is formed when the upstream plasma is slowed by the backpressure of the accelerated particles and the overall compression ratio changes from the test-particle value as a result of high energy particles escaping and/or a change in the shocked plasma’s adiabatic index $`\mathrm{\Gamma }`$. As indicated by our $`\gamma _0=1.4`$ example (Section 4.1), mildly relativistic shocks act as nonrelativistic ones showing a dramatic weakening of the subshock combined with a large increase in $`r`$ (Fig. 8). These changes result in a particle distribution which is both steeper than the test-particle power law at low momenta and flatter at high momenta (Fig. 9).

In faster shocks (i.e., $`\gamma _03`$), the initial test-particle spectrum is steep enough that particle escape is unimportant so only changes in $`\mathrm{\Gamma }`$ determine $`r`$ (Eq. 16). In contrast to nonrelativistic shocks where the production of relativistic particles causes the compression ratio to increase, we show that $`r`$ decreases smoothly to $`3`$ as $`\gamma _0`$ increases and the fraction of fully relativistic shocked particles approaches one (Fig. 4).

Our most important result is that efficient, mildly relativistic shocks do not produce particle spectra close to the so-called ‘universal’ power law having $`\sigma 4.3`$. This may be important for the interpretation of gamma-ray bursts since the internal shocks assumed responsible for converting the bulk kinetic energy of the fireball into internal particle energy may be mildly relativistic and the external shocks, believed responsible for producing gamma-ray burst afterglows, will always go through a mildly relativistic phase (see Piran, 1999, for a comprehensive review of gamma-ray bursts).

The authors thank M. G. Baring, L. O’C. Drury, F. C. Jones, and E. Parizot for helpful discussions. D.C.E. is grateful to the Department of Physics and Astronomy at Rice University and to the Dublin Institute for Advanced Studies for hosting visits where part of this work was done.


text:
# Observation of room-temperature spontaneous phase segregation in overdoped Bi2Sr2CaCu2O8+x single crystals

Correspondence author, \]email: lilu@aphy.iphy.ac.cn.

## Abstract

The occurrence and development of phase inhomogeneity in Bi2Sr2CaCu2O8+x single crystals have been investigated by measuring their superconducting transition patterns using a high-sensitivity ac magnetometer. We find that the overdoped Bi2Sr2CaCu2O8+x single crystals, even if they were single-phased immediately after a high-temperature annealing process, show progressively pronounced multiple superconducting transitions with time as they were kept at room temperature in a dry air atmosphere. The results reveal that the oxygen dopants in the crystal are still mobile at room temperature. These dopants tend to re-arrange and thus form different phases, at a characteristic time scale of one or two weeks or so. Moreover, the diamagnetization loops in the newly segregated phase are very weak, presumably due to the existence of additional inhomogeneity down to a nanometer scale. Our results also show that the optimally doped Bi2Sr2CaCu2O8+x crystals seem capable to remain in a single phase with time.

preprint: submitted to Phys. Rev. B

That the cuprates are of a single phase at a given temperature and a carrier concentration is an essential assumption in the widely-used phase diagrams, phasediagrams1 ; phasediagrams2  based on which many experimental aspects in the normal and superconducting states have been correlated and discussed. However, it is often a primary difficulty to obtain single crystals with high structural perfection and doping uniformity, yet with large enough size suitable for many experiments. Besides the electronic phase separation,stripe1 ; stripe2  chemical phases separation at mesoscopic or even larger scales, often related to the inhomogeneous distribution or domain formation of oxygen vacancies or dopants, has also been observed in many families of the cuprates,review1 ; review2  for examples, in La2CuO4+x, LCO4x1 ; LCO4x2 ; LCO4x3 ; LCO4x4  in lightly Sr-doped La2-xSrxCuO4,LSCOcho  in Nd2-xCexCuO4-y, Skelton  and in underdoped YBa2Cu3o7-x where vacancy ordering domains of different types form at the chain layers.YBCO1 ; Chen\_review  Phase inhomogeneity has also been observed in Bi2Sr2CaCu2O8+x by low temperature scanning tunneling microscopic studies. Pan ; Lang  With all these observations, however, one still cannot rule out the possibility that the phase separation is an extrinsic effect, avoidable if the crystals are grown up and annealed properly. In this paper, through investigating the time evolution of the diamagnetization signal, we show that the oxygen dopants in originally single-phased overdoped Bi2Sr2CaCu2O8+x single crystals tend to redistribute spontaneously at room temperature, which eventually and unavoidably leads to the occurrence of phase separation and inhomogeneity in the titled compound.

The Bi2Sr2CaCu2O8+x single crystals used in this experiment were proven to be of high structural quality and stoichiometric purity with the resolution limits of our X-ray diffraction and energy dispersion X-ray studies. To further guarantee the uniformity of the crystal, only those tiny pieces of crystalline sheets (typically 0.5 mm $`\times `$ 0.5 mm $`\times `$ 0.05 mm, $``$ 10-2 mm3 in volume) cleaved and cut from the center of larger uniform ones were used for investigation. Each crystal was carefully annealed to maintain a certain and uniform level of overdoping. Then, its susceptibility was measured in a weak ac magnetic field (0.3 to 60 G) using a high-sensitivity vibration-sample magnetometer specially designed for measuring sub-millimeter-sized specimens.setup  This magnetometer uses a double-synchronous detection technique, with a high-frequency ac magnetic field at 10 kHz (at which the Bi2Sr2CaCu2O8+x crystal is fully penetrated in the normal state), and a low-frequency vibration of the sample at 4.4 Hz. It reaches a resolution of $`1\times 10^8`$ emu in susceptibility (or $`3\times 10^{10}`$ emu in magnetic moments) for a sample of volume 1 mm3 in a field of 30 G. With such a sensitivity we can resolve any minority superconducting phase whose total volume is as small as 10-7 mm3, at the same time without suppressing the superconductivity. Our measurement on optimally doped Bi2Sr2CaCu2O8+x crystals reveals a single sharp transition (will be shown in Fig. 5), which proves the reliability of this home-made setup on the one hand and the high quality of the pristine crystals used in this experiment on the other hand.

Figure 1 shows the superconducting transitions of two overdoped Bi2Sr2CaCu2O8+x single crystals, samples #1 and #1’, measured immediately after an annealing treatment in O2 atmosphere at 550 oC for 52 hours. Sample #1 is measured in the configuration of B$``$c, whereas sample #1’ is a small piece cut from sample #1, in order to be mounted into the quartz-tube sample holder in the B$``$ab configuration. The 10% to 90% transition width for sample #1’ is $``$ 1 K in a field of 0.6 G perpendicular to the ab-plane, which demonstrates that an overdoped crystal could be in a uniform single phase.

For the case of B$``$c, the diamagnetization signal at low temperatures reflects the nature of intrinsic Josephson couplings along the c-direction. cJosephson  To distinguish the transition regime from this Josephson coupling regime, we find it helpful to plot the data in a logarithmic scale (the inset of Fig. 1). In such a plot the temperature dependence of the susceptibility is very steep in the transition regime but relatively “flat” in the Josephson coupling regime. The crossover between the two regimes takes place at $``$ 82.5 K for sample #1.

Although the overdoped crystals show a single superconducting transition if measured immediately after an annealing treatment, multiple transitions at different levels are always seen after the crystals have been kept at room temperature for a period of time. In our experiment, nine overdoped crystals have been examined, seven of them eventually exhibited a multiple transition pattern, and two of them underwent significant broadening in transition width, after they have been stored for a certain period of time.

Figure 2 shows the time evolution of the transition pattern for sample #2. This crystal was annealed in O2 atmosphere at 450 oC for 48 hours. It was afterwards kept at room temperature in a dry air environment as protected by silica gel, which is a common way of sample storage in most of the laboratories. In a three-month period the susceptibility of the crystal in the B$``$c configuration was measured for several times. From the panel (b) of Fig.2 it can be seen that the 88 K phase grew up significantly with time, and reached $``$ 10% in diamagnetization at the end of the three-month period.

The measurement performed in the B$``$ab configuration on sample #2’ (which is a small piece taken from sample #2 at the end of the three-month period) also revealed the formation of a very pronounced secondary phase, as shown in the panel (a) of Fig. 2.

The growing up in volume of a higher-Tc phase should reflect the change in local oxygen content and/or arrangement in the crystals. Only with such a mechanism can then the change be reversible, i.e., a single sharp transition can be restored by high-temperature annealing — which was exactly the case for our crystals. Our results therefore reveal that the balanced oxygen content in overdoped Bi2Sr2CaCu2O8+x is a function of the temperature and tbe oxygen partial pressure of the crystal’s environment. The nearly balanced oxygen content in a crystal, as maintained by high-temperature annealing in O2-atmosphere, must become over saturated as the crystal is cooled down to room temperature. Therefore, the oxygen dopants tend to redistribute and even to escape from the crystal, which finally creates the higher-Tc phase. For this reason, the newly formed phase is likely near the edge of the crystals. It should occupy a substantially large volume, in order to cause an apparent diamagnetic signal of $``$10% for sample #2, and up to $``$ 50% for sample #2’.

The distinct multiple-transition pattern indicates the formation of two types of domains in the ab-plane of the crystal, each having a different oxygen content. The mechanism of forming such domains might be similar to that of the vacancy ordering domains in underdoped YBa2Cu3O7-δ.YBCO1 ; Chen\_review  It has been shown through numerical simulations numerical  that an ordered distribution of oxygen vacancies commensurable with the lattice, i.e., forming some types of staging, is energetically more favored than a random distribution.

A non-trivial mobility for the oxygen dopants at room temperature is another key factor for the occurrence of spontaneous phase separation. It has been shown in La2CuO4+x that the oxygen mobility is linked to the lattice imperfections, e.g., it is strongly enhanced by the presence of planar defects.mobility  For Bi2Sr2CaCu2O8+x samples, the existence of incommensurate modulation along the b-axis incommensurate1 ; incommensurate2  might in some way influence the oxygen mobility, which needs to be further studied.

Figure 3 shows the field dependence of the transition pattern of sample #2’ three months after the annealing, and sample #2 24 days after the annealing. Two transitions can be recognized in weak magnetic fields, one onsets at $``$ 81.6 K and the other above 86 K. The nucleating and growing up of the superconducting domains below 86 K build up substantially large or enormous amount of diamagnetization loops in the ab-plane of the crystal. However, these loops are very weak, as they can be drastically suppressed by a field of $``$ 1 Gauss, and almost washed out by a field of a few tens Gauss. Such a sensitive field dependence should come from a Josephson-like weak superconductivity, reflecting the existence of additional inhomogeneity in the ab-plane of the newly segregated phase down to the coherence length scale. Other\_field\_dependence

The diffusion process of the dopants to a lower-energy distribution is presumably driven by the thermal energy which is significantly lower at room temperature compared to that in high temperature annealing. This factor makes the diffusion process hardly completed, givn rise to the additional inhomogeneity in the ab-plane. This picture seems consistent with the recent observation of energy gap inhomogeneity at nanometer scale using scanning tunneling microscopy.Pan ; Lang  With the best single crystals they obtained, Pan and co-workers found that the superconducting energy gap of Bi2Sr2CaCu2O8+x compound is inhomogenous from place to place down to the nanometer scale, though the lattice is highly perfect.

Finally, Fig. 5 shows the transition patterns of two optimally doped Bi2Sr2CaCu2O8+x crystals (the 88 K phase) 15 days (sample #3) and six months (sample #4) after being annealed in Ar atmosphere (600 oC, two weeks). These crystals are from the same batch as sample #1 and #2. It seems that the tendency of spontaneous phase separation at room temperature is much insignificant for optimally doped crystals, as their transition pattern remains nearly the same after the six months period, and shows no significant field dependence at a level of a few Gauss. It indicates that the oxygen dopants in optimally doped crystals are near to a minimum energy distribution which is homogenous down to the coherence length scale.

###### Acknowledgements.

One of us (LL) thanks H. H. Wen and S. H. Pan for helpful discussion. This work was supported by the National Science Foundation of China, and by the National Center for R & D of Superconductivity, China.