text:
# A New Interpretation of Flux Quantization

## I Introduction

Flux quantization is one of the most fascinating properties of superconducting states. This phenomenon was found by Deaver and Fairbank1 and Doll and Näbauer2 in 1961. They observed that the magnetic flux trapped in a superconducting cylinder is quantized in units of hc/2e. Recently, half-integer flux quantization has attracted much attention in connection with the pairing symmetry in high Tc cuprates.3

Byers and Yang4 proposed that flux quantization follows because the free energy of the superconducting state is periodic in the unit of the flux, hc/2e, if electrons are paired. In fact, they assumed that i) the free energy, $`F_N`$, of the normal state is independent of the flux, i.e., $`F_N(\mathrm{\Phi })=\mathrm{constant}`$, and ii) the pairing energy is constant as a function of the flux.5 Consequently, they ascribed flux quantization to the minimum of the kinetic energy of the superconducting state at $`\mathrm{\Phi }/(hc/2e)=\mathrm{integer}`$. Based on the theory of Byers and Yang, Little and Parks6 pointed out that the transition temperature $`T_c`$ is also a periodic function of the enclosed flux $`\mathrm{\Phi }`$, which is called Little-Parks oscillations.

However, recent discovery of persistent current in normal metal rings at low temperature casts a doubt on the first assumption, that the free energy of the normal state is essentially independent of the flux.7-9 Since the electronic states of normal metal rings are periodic in flux $`\mathrm{\Phi }`$, with period hc/e, the persistent current is given by $`I(\mathrm{\Phi })=c\times F_N(\mathrm{\Phi })/\mathrm{\Phi }`$. Therefore, the first assumption is not compatible with the experimental observation of the persistent current in normal metal rings.

In this paper we report that flux quantization and the Little-Parks oscillations are caused by the flux dependence of the pairing energy. In addition, we show that the flux dependence of the kinetic energy does not lead to any significant oscillations in the transition temperature $`T_c`$ and the energy gap parameter $`\mathrm{\Delta }_0`$ for superconducting cylinders. As the flux increases, the pairing matrix element is reduced since we must pair single-particle states which are not degenerate and have different density distributions.10,11 When the flux $`(\mathrm{\Phi }`$) is equal to $`\mathrm{\Phi }=n\times hc/2e`$ (with integer n), the pairing interaction and the free energy of the superconductor remain unchanged (even n) or almost unchanged (odd n) from their values with zero flux, since we can pair degenerate partners, with the same density distribution, due to the energy level crossing. Accordingly, we find flux quantization and Little-Parks oscillations. It is interesting to note that Schwartz and Cooper12 obtained a similar result for generalized pairing with different angular momentum in superconducting cylinders.

## II Theory of Byers and Yang

To study the effect of the Aharonov-Bohm flux $`\mathrm{\Phi }`$ on the superconducting states in cylinders, Byers and Yang4 considered the variation of the energy levels E of a single particle in the presence of the flux:

$$E=\frac{1}{2m}[p_r^2+p_z^2+\frac{\mathrm{}^2}{r^2}(m+\frac{e}{hc}\mathrm{\Phi })^2].$$
(1)
Here $`p_r`$ and $`p_z`$ denote momenta in the radial and z directions. Employing the BCS pairing between states m and -m, they noted that the average (kinetic) energy for the states m and -m, near $`\mathrm{\Phi }=+0`$, increases with $`\mathrm{\Phi }`$ like

$$\mathrm{constant}+\frac{\mathrm{}^2}{2\mathrm{m}\mathrm{r}^2}(\frac{\mathrm{e}\mathrm{\Phi }}{\mathrm{hc}})^2.$$
(2)
Furthermore, at $`2e\mathrm{\Phi }/hc=1`$, since pairing between

$$m+\frac{e}{hc}\mathrm{\Phi }=m+\frac{1}{2}=\frac{1}{2}\mathrm{and}\frac{1}{2},\frac{3}{2}\mathrm{and}\frac{3}{2},\mathrm{etc}.,$$
(3)
is preferred, the kinetic energy per particle remains the same as for the case $`\mathrm{\Phi }=0`$. Then the additional energy for each of these pairs, in the neighborhood of $`2e\mathrm{\Phi }/hc`$, increases with $`\mathrm{\Phi }`$ in a parabolic manner:

$$2\times (\frac{\mathrm{}^2}{2mr^2})[(\frac{e\mathrm{\Phi }}{hc})\frac{1}{2}]^2.$$
(4)
Consequently, the kinetic energy of the Cooper pairs has the minimum at $`\mathrm{\Phi }=n\times hc/2e`$ (with integer n). Based upon this observation, Byers and Yang4 asserted that flux quantization is due to the flux-dependent minimum of the kinetic energy of the Cooper pairs. On the other hand, the corresponding BCS wavefunctions at $`\mathrm{\Phi }=n\times hc/2e`$ are given by13

$$(u_\nu +v_\nu a_{\nu +m}^+a_{\nu +m}^+)|0>;(n\mathrm{even})$$
(5)
$$(u_{\nu +1/2}+v_{\nu +1/2}a_{\nu +m+1}^+a_{\nu +m}^+)|0>;(n\mathrm{odd})$$
(6)
with $`n=2m(n=2m+1)`$ for even (odd) values of n. The operator $`a_{\nu +m}^+`$ creates an electron in the state $`(\nu +m)`$.

However, the kinetic energy of the normal electrons will show the same flux dependence, as confirmed by the observation of the persistent current in normal metal rings.7-9 Therefore, we expect that this flux dependence will be canceled in the superconducting state. More precisely, using the BCS theory, we show that the variation of the kinetic energy with the flux does not lead to any oscillations of the condensation energy (i.e., $`\mathrm{\Delta }_0`$) and the free energy (i.e., $`T_c`$) of superconducting cylinders.

Assuming a constant pairing matrix element V,4,5 one finds the BCS gap equation at $`T=0K`$ and the $`T_c`$ equation:14-16

$$\mathrm{\Delta }_m=V\underset{m^{}}{}\frac{\mathrm{\Delta }_m^{}}{2\sqrt{(\frac{ϵ_{m^{}\mu }+ϵ_{m^{}\mu }}{2})^2+\mathrm{\Delta }_m^{}^2}}$$
(7)
$$1=V\underset{m}{}\frac{tanh\frac{ϵ_{m\mu }}{2T_c}+tanh\frac{ϵ_{m\mu }}{2T_c}}{2(ϵ_{m\mu }+ϵ_{m\mu })}$$
(8)
where $`ϵ_{m\mu }`$ is the one-particle energy (measured from the Fermi energy $`E_F`$) with $`\mu =e\mathrm{\Phi }/hc`$. Contrary to Byers and Yang, it is clear that the kinetic energy does not give rise to any significant oscillations in the gap parameter, $`\mathrm{\Delta }_m=\mathrm{\Delta }_0`$, since its effect is only to shift the BCS cutoff region from $`\omega _D<ϵ<\omega _D`$ to $`\omega _D+ϵ_\mu <ϵ<\omega _D+ϵ_\mu `$ (where $`ϵ_\mu \frac{\mathrm{}^2\mu ^2}{2mr^2}`$):

$$\mathrm{\Delta }_0=V\underset{m}{}\frac{\mathrm{\Delta }_0}{2\sqrt{(ϵ_m+ϵ_\mu )^2+\mathrm{\Delta }_0^2}}.$$
(9)
Use of $`ϵ_{m\mu }+ϵ_{m\mu }=2(ϵ_m+ϵ_\mu )`$ has been made. For the $`T_c`$, Eq. (8) shows indeed a minor decrease in $`T_c`$ (for nanoscale systems) due to the flux, which is inconsistent with the ($`T=0K`$) gap equation, Eq. (7), and has been elusive.5,6 However, this decrease is completely negligible for bulk systems and it originates from the exclusion of the Cooper pairs with negative excitation energy.17,18 Since analytic solution is not available, we rely on the numerical calculations and merely emphasize that the correction is in the denominator:

$`1`$ $`=`$ $`V{\displaystyle \underset{m}{}}{\displaystyle \frac{tanh\frac{ϵ_{m\mu }}{2T_c}+tanh\frac{ϵ_{m\mu }}{2T_c}}{2(ϵ_m+ϵ_\mu )}}`$ (10)
$`=`$ $`V{\displaystyle \underset{m}{}}{\displaystyle \frac{1}{2(ϵ_m+ϵ_\mu )}}{\displaystyle \frac{sinh\frac{ϵ_m+ϵ_\mu }{T_c}}{cosh\frac{ϵ_{m\mu }}{2T_c}cosh\frac{ϵ_{m\mu }}{2T_c}}}.`$ (11)
Note that for bulk systems m is of the order of 1,000 near $`E_F`$, while $`\mu `$ is of the order of 1.

For example, we solve Eqs. (7) and (8) numerically for a superconducting cylinder with the inner radius $`r_a=7,200\AA `$, the outer radius $`r_b=8,000\AA `$, and the height $`L=500\AA `$. For a thin cylinder, i.e. $`r_br_a<<r_b`$, the electron energy is approximately written as $`ϵ_{m\mu }\frac{\mathrm{}^2}{2m}[\frac{(m\mu )^2}{r_ar_b}+(\frac{\mathrm{}\pi }{r_br_a})^2+(\frac{n_z\pi }{L})^2]`$.19,20 (More calculational details will be explained in the next section.) We assumed $`E_F=1.0eV`$ and $`\omega _D=0.02eV`$. There are 1,286,200 states in the BCS cutoff range. Figure 1 shows $`\mathrm{\Delta }_0`$ and $`T_c`$ versus $`\mu =e\mathrm{\Phi }/hc`$. As can be seen, there is no change in $`\mathrm{\Delta }_0`$ (i.e., the condensation energy) and $`T_c`$ (i.e., the free energy of superconducting state) due to the variation of the kinetic energy with the flux $`\mathrm{\Phi }`$. For a certain V ($`=17.0/\mathrm{volume}`$ \[eV\]), the gap parameters $`\mathrm{\Delta }_0`$ are the same 2.144 886 meV for $`\mu =0`$ and $`\mu =0.25`$, whereas the $`T_c^{}s`$ are 1.212 888 meV ($`\mu =0`$) and 1.211 921 meV ($`\mu =0.25`$), respectively. It is noteworthy that the decrease in $`T_c`$ scales approximately $`1/\sqrt{N}`$, with N= number of states in the BCS cutoff range. Consequently, the $`T_c`$ decrease due to the exclusion of the Cooper pairs with negative excitation energy17,18 can be completely ignored in the bulk limit. Note also that $`2\mathrm{\Delta }_0/T_c3.537`$, in agreement with the bulk result.

## III Flux-dependent pairing interaction

Now we consider flux-dependence of the pairing matrix element and the pairing interaction in a superconducting cylinder whose axis coincides with the z-axis. The inner and outer radii are $`r_a`$ and $`r_b`$, and the height is L. First, we determine one-particle eigenfunctions and eigenvalues in the presence of the Aharonov-Bohm flux $`\mathrm{\Phi }`$. The Hamiltonian is given by

$$H_0=\frac{\mathrm{}^2}{2m_e}[\frac{^2}{r^2}+\frac{1}{r}\frac{}{r}\frac{1}{r^2}(i\frac{}{\varphi }\mu )^2+\frac{^2}{z^2}]$$
(12)
where $`\mu =e\mathrm{\Phi }/hc`$. Assuming a solution of the form

$$\varphi (r,\varphi ,z)=R(r)sin(\frac{n_z\pi z}{L})e^{im\varphi }.$$
(13)
we obtain the Bessel equation for $`R(r)`$:

$$R^{\prime \prime }+\frac{R^{}}{r}+[k^2\frac{(m\mu )^2}{r^2}]R=0.$$
(14)
From the boundary condition $`R(r=r_a)=R(r=r_b)=0`$, we find the radial solutions and eigenvalues:19,20

$$R_{m\mu }(r)=𝒩[N_{m\mu }(kr_a)J_{m\mu }(kr)J_{m\mu }(kr_a)N_{m\mu }(kr)]$$
(15)
$$E_{lmn_z}=\frac{\mathrm{}^2}{2m_e}[k_{ml}^2+(\frac{n_z\pi }{L})^2](n_z=1,2,3,\mathrm{},m=0,\pm 1,\pm 2,\mathrm{})$$
(16)
where $`𝒩`$ is the normalization constant and $`k_{ml}`$ is the $`l`$th root of the transcendental equation

$$N_{m\mu }(kr_a)J_{m\mu }(kr_b)J_{m\mu }(kr_a)N_{m\mu }(kr_b)=0.$$
(17)
From Eq. (16) we calculate the single-particle energy as a function of the flux $`\mathrm{\Phi }`$, shown in Figure 2, near $`E_F=1.0`$ eV for $`r_a=75\AA `$ and $`r_b=150\AA `$. Ignoring the z coordinate, we find 15 basis pair states within the BCS cutoff range for $`\omega _D=0.02eV`$. We consider a small size ring for simplicity and because of the difficulty in calculating the high order Bessel functions. As $`\mu `$ increases, some states near 1.02 eV get out of the BCS cutoff range, while some states below 0.98 eV are coming in the range. We disregard this complication, since this effect will be negligible in the bulk limit. Moreover, it is caused by the simple approximation, which can be avoided in the strong-coupling theory.21 Near $`\mu =0.5`$ we chose $`E_F0.99`$ eV, which is obtained from the $`E_F`$ for $`\mu =0`$ by applying the flux. Although the spectra near $`\mu =0`$ and $`\mu =0.5`$ are similar, they are not exactly the same due to the finite size effect. Notice the crossings of energy levels at $`\mathrm{\Phi }=hc/2e\times n`$ (with integer n).

Second, we derive the pairing matrix element for an Einstein phonon model with the phonon Green’s function $`D(xx^{})`$. Using the equivalent electron-electron potential in the electron-phonon problem,11,22,23

$$V(xx^{})V_0D(xx^{}),$$
(18)
with $`x=(𝐫,t)`$, we get the pairing matrix element:

$$V_{m,m^{}}=V_0_{r_a}^{r_b}R_{m^{}\mu }^{}(r)R_{m^{}\mu }^{}(r)R_{m\mu }(r)R_{m\mu }(r)r𝑑r;\frac{1}{2}+\delta <\mu <\frac{1}{2}\delta $$
(19)
$$V_{m,m^{}}=V_0_{r_a}^{r_b}R_{m^{}\mu }^{}(r)R_{m^{}+1\mu }^{}(r)R_{m+1\mu }(r)R_{m\mu }(r)r𝑑r;\frac{1}{2}\delta <\mu <1\delta $$
(20)
$$V_{m,m^{}}=V_0_{r_a}^{r_b}R_{m^{}+1\mu }^{}(r)R_{m^{}+1\mu }^{}(r)R_{m+1\mu }(r)R_{m+1\mu }(r)r𝑑r;1\delta <\mu <\frac{3}{2}\delta $$
(21)
where $`\delta 1/4`$. We have made use of the BCS wavefunction with the flux:

$$\stackrel{~}{\varphi }_{BCS}=\underset{m}{}(u_m+v_ma_{m\mu }^+a_{m\mu }^+|0>;\frac{1}{2}+\delta <\mu <\frac{1}{2}\delta $$
(22)
$$\stackrel{~}{\varphi }_{BCS}=\underset{m}{}(u_m+v_ma_{m\mu }^+a_{m+1\mu }^+|0>;\frac{1}{2}\delta <\mu <1\delta $$
(23)
$$\stackrel{~}{\varphi }_{BCS}=\underset{m}{}(u_m+v_ma_{m\mu }^+a_{m+2\mu }^+|0>;1\delta <\mu <\frac{3}{2}\delta .$$
(24)
Finally, we obtain the BCS gap equation14-16

$$\mathrm{\Delta }_m=\underset{m^{}}{}V_{mm^{}}\frac{\mathrm{\Delta }_m^{}}{4\sqrt{(\frac{ϵ_{m^{}\mu }+ϵ_{m^{}\mu }}{2})^2+\mathrm{\Delta }_m^{}^2}}(tanh\frac{E_{m^{}\mu }}{2T}+tanh\frac{E_{m^{}\mu }}{2T})$$
(25)
where the excitation energy $`E_{m\mu }`$ is

$$E_{m\mu }=\frac{1}{2}(ϵ_{m\mu }ϵ_{m\mu })+\sqrt{(\frac{ϵ_{m\mu }+ϵ_{m\mu }}{2})^2+\mathrm{\Delta }_m^2}.$$
(26)
Since the analytic expression of the matrix element is not available, we rely on the numerical method to calculate $`T_c`$ from Eq. (24). Note that Eq. (24) leads to the negative excitation energy, when $`ϵ_{m\mu }`$ and $`ϵ_{m\mu }`$ have different signs near $`E_F`$, as in the cases of a superconductor with a uniform exchange field17 and the current-carrying superconducting state.18 Nevertheless, in this case the possibility of the negative excitation energy is not significant as shown in Fig. 1, because $`m1,000`$ near $`E_F`$ while $`\mu 1`$ in the bulk limit. In other words, the exclusion of the contribution of the negative excitation does not lead to any noticeable effect in $`T_c`$, as shown in Fig. 1.

## IV Explanation of Flux Quantization

Flux quantization may be explained from the flux-dependence of the pairing matrix element $`V_{mm^{}}`$. We stress the periodicity of the matrix element with the flux:

$$V_{mm^{}}(\mathrm{\Phi }=\frac{hc}{2e}\times \mathrm{even}\mathrm{n})=\mathrm{V}_{\mathrm{mm}^{}}(\mathrm{\Phi }=0).$$
(27)
$$V_{mm^{}}(\mathrm{\Phi }=\frac{hc}{2e}\times \mathrm{odd}\mathrm{n})\mathrm{V}_{\mathrm{mm}^{}}(\mathrm{\Phi }=0).$$
(28)
This behavior is caused by the energy level crossings at $`\mathrm{\Phi }=hc/2e\times n`$ (or $`\mu =\mathrm{integer}\mathrm{or}\mathrm{half}\mathrm{integer}`$). It was shown that the level crossings (at integer $`\mu `$) are due to the cylindrical symmetry of the system.4,20 Since for odd n the energy level crossing occurs at lower energy than for even n, the matrix elements are slightly different for even and odd n’s.

In the neighborhood of $`\mathrm{\Phi }=0+`$ we pair states ($`m\mu `$) and ($`m\mu `$) which are not degenerate and have a different density distribution.10,11 Accordingly, the pairing matrix element $`V_{mm^{}}`$ is decreasing with the flux, compared to $`V_{mm^{}}`$ at zero flux. For a large cylinder, we may employ the asymptotic radial wavefunction:24

$$R_{m\mu }(r)=𝒩\{[1+a_{m\mu }\frac{1}{r_ar}b_{m\mu }(\frac{1}{2r^2}+\frac{1}{2r_a^2})]\frac{sink_l(rr_a)}{\sqrt{r}}\}$$
(29)
where

$$a_{m\mu }=[\frac{4(m\mu )^21}{8k_l}]^2,b_{m\mu }=\frac{[4(m\mu )^21][4(m\mu )^29]}{(8k_l)^2},$$
(30)
and

$$k_l=\frac{l\pi }{r_br_a}.$$
(31)
It is obvious that the algebraic terms lead to the different density distributions of states $`(m\mu )`$ and $`(m\mu )`$ and the decrease of the matrix element:

$$V_{mm^{}}(\mathrm{\Phi })[1C^{}(\frac{2m^3}{k_l^2})^2\mu ^2][1C^{}(\frac{2m^3}{k_l^{}^2})^2\mu ^2]\times V_{mm^{}}(\mathrm{\Phi }=0),$$
(32)
where $`C^{}(1/r_ar_b1/2r_b^21/2r_a^2)^2`$.

Figure 3 shows the radial density distributions $`|R(r)|^2r`$ of the $`m=\pm 54`$ states when $`\mu =0.25`$, (i.e., $`m\mu =53.75,m\mu =54.25`$) for the same ring as Fig. 2. The solid line corresponds to $`m\mu =54.75`$, whereas the dotted line corresponds to $`m\mu =54.25`$, in which the electron is rotating faster, due to bigger angular momentum projection along the z-axis, leading to the shift of the radial distribution towards the outer radius. The difference in the radial density distributions then induces the decrease of the pairing matrix element with increasing $`\mu `$. Figure 4 shows the matrix element $`V_{mm^{}}V_{ij}`$ for $`\mu =0`$. Notice that it is $`15\times 15`$ symmetric matrix, since there are 15 pair states in the BCS cutoff range. The decrease of the matrix element with $`\mu `$ is not big enough to be noticeable in Fig. 4, though. We just mention the decrease of biggest matrix element $`V_{m=53,m^{}=51}V_{i=2,j=3}`$, between pair states with the same $`l`$ value, from 0.0002450 to 0.0002448.

This observation and the periodicity of $`V_{mm^{}}V_{i,j}`$ prove that flux quantization is due to the minimum of the pairing energy of the superconducting cylinders as a function of the flux. Accordingly, the free energy of the superconducting state also has minimum at $`\mathrm{\Phi }=hc/2e\times n`$ since the contribution of the kinetic energy is constant. In other words, pairing between the degenerate states caused by the energy level crossing at $`\mathrm{\Phi }=hc/2e\times n`$ is the physical origin of flux quantization. Thus the Little-Parks oscillations follow accordingly. From Eqs. (24) and (31) we may estimate roughly the decreases of the energy gap at $`T=0K`$ and $`T_c`$:

$`\mathrm{\Delta }_m(\mathrm{\Phi })`$ $``$ $`\mathrm{\Delta }_m(\mathrm{\Phi }=0)(18C^{}<{\displaystyle \frac{m^6}{k_l^4}}>\mu ^2/N_0V_0),`$ (33)
$`T_c`$ $``$ $`T_{c0}(18C^{}<{\displaystyle \frac{m^6}{k_l^4}}>\mu ^2/N_0V_0),`$ (34)
where $`N_0`$ is the density of states at $`E_F`$ and $`<>`$ denotes the average over the states in the BCS cutoff range. $`T_{c0}`$ is the transition temperature in the absence of the magnetic flux.

Figure 5 shows the $`T=0K`$ anisotropic gap parameter $`\mathrm{\Delta }_i`$ for $`\mu =0`$, calculated from Eq. (24) by employing the matrix element $`V_{i,j}V_{mm^{}}`$ as in Fig. 4. Even though the coupling between the second and third states is strongest, we find the biggest gap at the tenth state, i.e., $`\mathrm{\Delta }_{i=10}`$, due to its strong couplings to several neighbor states.

Figure 6 displays the decrease of the (average) energy gap with $`\mu `$. The solid line denotes the ($`T=0K`$) average gap parameter, $`\overline{\mathrm{\Delta }}=_i\mathrm{\Delta }_i/15`$. It should be noticed that this decrease, so much enhanced as a result of the finite size effect, persists even in the bulk limit, as Eq. (32) shows. The different variations near $`\mu =0`$ and $`\mu =0.5`$ are also due to the finite size effect. The dashed line shows the isotropic gap parameter, according to the assumptions of Byers and Yang.4 It exhibits a small decrease with the flux too, which quickly disappears with the increase of the system size. For the isotropic solution, we included $`V_{ii}=V`$ besides $`V_{ij}=V=`$ constant. Therefore, we get the bigger gap. The semicircle is the $`T_c`$ when $`\mu =0`$. Since, as the flux is increased, the finite temperature gap equation, Eq. (24) eliminates the contributions from the Cooper pairs with negative excitation energy,17,18 $`T_c`$ drops to zero immediately. If we consider the Cooper problem,25 it is easy to understand why we need to eliminate the Cooper pair basis states with different signs in energy near $`E_F`$. However, this behavior is greatly amplified because of the finite size effect and it will not cause any problem in realistic systems, as shown in Fig. 1. Furthermore, the strong coupling theory21 will reduce the effect even more. For nanoscale systems it may be desirable to use the canonical ensemble approach.26 We expect that for a reasonably large system, we find a continuous (small) variation of the $`T_c`$, consistent with the variation of the energy gap.

## V Discussion

It is clear that more study is needed to understand the details of flux quantization. In particular, it is essential to consider a big hollow cylinder for comparison with the experiment, which requires high order Bessel functions. This study indicates that the parabolic background5,6 in the $`T_c`$ oscillations may be due to the magnetic field penetration into the superconductor, which is not addressed here. It is also interesting that disordered $`Au_{0.7}In_{0.3}`$ cylinders27 show half-integer flux quantization as in the case of high $`T_c`$ cuprates.3

This study also implies that the exact calculation of the eigenenergies using high order Bessel functions and the magnetic field penetration into the metals are important in calculating the persistent current in normal metals.

Finally, flux quantization plays an important role in the vortex state28 which may be understood by generalizing this approach.

## VI Conclusion

We have shown that flux quantization and the Little-Parks oscillations in superconducting cylinders are due to the flux dependence of the pairing interaction, while the flux dependence of the kinetic energy is canceled as evidenced by the observation of the persistent current in normal metals. When the flux $`\mathrm{\Phi }`$ is given by $`\mathrm{\Phi }=hc/2e\times n`$ (with integer n), the pairing interaction and the free energy are unchanged (even n) or almost unchanged (odd n), due to the degenerate-state pairing, compared to those in the absence of the flux. Accordingly, flux quantization and the Little-Parks oscillations are obtained.

ACKNOWLEDGMENTS

M.P. thanks NSF-EPSCOR (Grant No. EPS9874782) for financial support. YJK is grateful to Faculty of Arts and Sciences at UPR-Mayaguez for release time.


format: text

text:
# The new analysis of the KLOE data on the ϕ→𝜂𝜋⁰𝛾 decay

## I Introduction

The lightest scalar mesons $`a_0(980)`$ and $`f_0(980)`$, discovered more than thirty years ago, became the hard problem for the naive quark-antiquark ($`q\overline{q}`$) model from the outset. Really, on the one hand the almost exact degeneration of the masses of the isovector $`a_0(980)`$ and isoscalar $`f_0(980)`$ states revealed seemingly the structure similar to the structure of the vector $`\rho `$ and $`\omega `$ mesons , and on the other hand the strong coupling of $`f_0(980)`$ with the $`K\overline{K}`$ channel pointed unambiguously to a considerable part of the strange quark pair $`s\overline{s}`$ in the wave function of $`f_0(980)`$. It was noted late in the 1970s that in the MIT bag model there are light four-quark scalar states and suggested that $`a_0(980)`$ and $`f_0(980)`$ might be these states . From that time $`a_0(980)`$ and $`f_0(980)`$ resonances came into beloved children of the light quark spectroscopy, see, for example, Refs. .

Ten years later there was proposed in Ref. to study radiative $`\varphi `$ decays $`\varphi a_0\gamma \eta \pi ^0\gamma `$ and $`\varphi f_0\gamma \pi ^0\pi ^0\gamma `$ to solve the puzzle of the lightest scalar mesons. Over the next ten years before the experiments (1998), this question was examined from different points of view .

Now these decays have been studied not only theoretically but also experimentally. The first measurements have been reported by the SND  and CMD-2  Collaborations which obtain the following branching ratios

$$Br(\varphi \gamma \pi ^0\eta )=(8.8\pm 1.4\pm 0.9)\times 10^5\text{[13]},$$

$$Br(\varphi \gamma \pi ^0\pi ^0)=(12.21\pm 0.98\pm 0.61)\times 10^5\text{[14]}$$

$$Br(\varphi \gamma \pi ^0\eta )=(9.0\pm 2.4\pm 1.0)\times 10^5,$$

$$Br(\varphi \gamma \pi ^0\pi ^0)=(9.2\pm 0.8\pm 0.6)\times 10^5\text{[15]}.$$

More recently the KLOE Collaboration has measured

$$Br(\varphi \gamma \pi ^0\eta )=(8.51\pm 0.51\pm 0.57)\times 10^5\text{in}\eta \gamma \gamma \text{[16]},$$

$$Br(\varphi \gamma \pi ^0\eta )=(7.96\pm 0.60\pm 0.40)\times 10^5\text{in}\eta \pi ^+\pi ^{}\pi ^0\text{[16]},$$

$$Br(\varphi \gamma \pi ^0\pi ^0)=(10.9\pm 0.3\pm 0.5)\times 10^5\text{[17]},$$

in agreement with the Novosibirsk data  but with a considerably smaller error.

In this work we present the new analysis of the recent KLOE data on the $`\varphi \eta \pi ^0\gamma `$ decay . In contradistinction to , we

1) treat the $`a_0`$ mass $`m_{a_0}`$ as a free parameter of the fit;

2) fit the phase $`\delta `$ of the interference between $`\varphi a_0\gamma \eta \pi ^0\gamma `$ (signal) and $`\varphi \rho ^0\pi ^0\eta \pi ^0\gamma `$ (background) reactions;

3) use new more precise experimental values of the input parameters.

All formulas for the $`\varphi (a_0\gamma +\rho ^0\pi ^0)\eta \pi ^0\gamma `$ reaction taking the background into account are shown in Sec.II. The results of the 4 different fits are presented in Sec.III. A brief summary is given in Sec.IV.

## II The formalism of the $`\varphi a_0\gamma \eta \pi ^0\gamma `$ and $`\varphi \rho ^0\pi ^0\eta \pi ^0\gamma `$ reactions

In Ref. was shown that the process $`\varphi a_0\gamma \eta \pi ^0\gamma `$ dominates in the $`\varphi \eta \pi ^0\gamma `$ decay (see also , where it was predicted in four-quark model). This was confirmed in . Nevertheless, the main background process $`\varphi \rho \pi ^0\eta \pi ^0\gamma `$ should be taken into account also (see ).

The amplitude of the background process $`\varphi (p)\pi ^0\rho ^0\gamma (q)\pi ^0(k_1)\eta (k_2)`$ is :

$$M_B=\frac{g_{\varphi \rho \pi }g_{\rho \eta \gamma }}{D_\rho (pk_1)}\varphi _\alpha k_{1\mu }p_\nu ϵ_\delta (pk_1)_\omega q_ϵϵ_{\alpha \beta \mu \nu }ϵ_{\beta \delta \omega ϵ}.$$
(1)
According to the one-loop mechanism of the decay $`\varphi K^+K^{}\gamma a_0`$, suggested in Ref., the amplitude of the signal $`\varphi \gamma a_0\gamma \pi ^0\eta `$ has the form:

$$M_a=g(m)\frac{g_{a_0K^+K^{}}g_{a_0\pi \eta }}{D_{a_0}(m)}\left((\varphi ϵ)\frac{(\varphi q)(ϵp)}{(pq)}\right),$$
(2)
where $`m^2=(k_1+k_2)^2`$, $`\varphi _\alpha `$ and $`ϵ_\mu `$ are the polarization vectors of $`\varphi `$ meson and photon, the forms of $`g_R(m)`$ and $`g(m)=g_R(m)/g_{RK^+K^{}}`$ everywhere over the $`m`$ region are in Refs.  and  respectively:

For $`m<2m_{K^+}`$

$`g(m)={\displaystyle \frac{e}{2(2\pi )^2}}g_{\varphi K^+K^{}}\{1+{\displaystyle \frac{1\rho ^2(m^2)}{\rho ^2(m_\varphi ^2)\rho ^2(m^2)}}\times `$ (3)
$`[2|\rho (m^2)|\mathrm{arctan}{\displaystyle \frac{1}{|\rho (m^2)|}}\rho (m_\varphi ^2)\lambda (m_\varphi ^2)+i\pi \rho (m_\varphi ^2)`$ (4)
$`(1\rho ^2(m_\varphi ^2))({\displaystyle \frac{1}{4}}(\pi +i\lambda (m_\varphi ^2))^2`$ (5)
$`\left(\mathrm{arctan}{\displaystyle \frac{1}{|\rho (m^2)|}}\right)^2)]\},`$ (6)
where

$$\rho (m^2)=\sqrt{1\frac{4m_{K^+}^2}{m^2}};\lambda (m^2)=\mathrm{ln}\frac{1+\rho (m^2)}{1\rho (m^2)};\frac{e^2}{4\pi }=\alpha =\frac{1}{137}.$$
(7)
For $`m2m_{K^+}`$

$`g(m)={\displaystyle \frac{e}{2(2\pi )^2}}g_{\varphi K^+K^{}}\{1+{\displaystyle \frac{1\rho ^2(m^2)}{\rho ^2(m_\varphi ^2)\rho ^2(m^2)}}\times `$ (8)
$`\times [\rho (m^2)(\lambda (m^2)i\pi )\rho (m_\varphi ^2)(\lambda (m_\varphi ^2)i\pi )`$ (9)
$`{\displaystyle \frac{1}{4}}(1\rho ^2(m_\varphi ^2))((\pi +i\lambda (m_\varphi ^2))^2(\pi +i\lambda (m^2))^2)]\}.`$ (10)
The mass spectrum is

$$\frac{d\mathrm{\Gamma }(\varphi \gamma \pi ^0\eta ,m)}{dm}=\frac{d\mathrm{\Gamma }_{a_0}(m)}{dm}+\frac{d\mathrm{\Gamma }_{back}(m)}{dm}+\frac{d\mathrm{\Gamma }_{int}(m)}{dm},$$
(11)
where the mass spectrum for the signal is

$$\frac{d\mathrm{\Gamma }_{a_0}(m)}{dm}=\frac{2}{\pi }\frac{m^2\mathrm{\Gamma }(\varphi \gamma a_0,m)\mathrm{\Gamma }(a_0\pi ^0\eta ,m)}{|D_{a_0}(m)|^2}=\frac{2|g(m)|^2p_{\eta \pi }(m_\varphi ^2m^2)}{3(4\pi )^3m_\varphi ^3}\left|\frac{g_{a_0K^+K^{}}g_{a_0\pi \eta }}{D_{a_0}(m)}\right|^2.$$
(12)
The mass spectrum for the background process $`\varphi \pi ^0\rho \gamma \pi ^0\eta `$ is :

$$\frac{d\mathrm{\Gamma }_{back}(m)}{dm}=\frac{(m_\varphi ^2m^2)p_{\pi \eta }}{128\pi ^3m_\varphi ^3}_1^1𝑑xA_{back}(m,x),$$
(13)
where

$`A_{back}(m,x)={\displaystyle \frac{1}{3}}{\displaystyle |M_B|^2}=`$ (14)
$`={\displaystyle \frac{1}{24}}(m_\eta ^4m_\pi ^4+2m^2m_\eta ^2m_\pi ^2\stackrel{~}{m_\rho }^22m_\eta ^4m_\pi ^2\stackrel{~}{m_\rho }^22m_\eta ^2m_\pi ^4\stackrel{~}{m_\rho }^2+`$ (15)
$`2m^4\stackrel{~}{m_\rho }^42m^2m_\eta ^2\stackrel{~}{m_\rho }^4+m_\eta ^4\stackrel{~}{m_\rho }^42m^2m_\pi ^2\stackrel{~}{m_\rho }^4+4m_\eta ^2m_\pi ^2\stackrel{~}{m_\rho }^4+m_\pi ^4\stackrel{~}{m_\rho }^4+`$ (16)
$`2m^2\stackrel{~}{m_\rho }^62m_\eta ^2\stackrel{~}{m_\rho }^62m_\pi ^2\stackrel{~}{m_\rho }^6+\stackrel{~}{m_\rho }^82m_\eta ^4m_\pi ^2m_\varphi ^22m^2m_\eta ^2m_\varphi ^2\stackrel{~}{m_\rho }^2+`$ (17)
$`2m_\eta ^2m_\pi ^2m_\varphi ^2\stackrel{~}{m_\rho }^22m^2m_\varphi ^2\stackrel{~}{m_\rho }^4+2m_\eta ^2m_\varphi ^2\stackrel{~}{m_\rho }^42m_\varphi ^2\stackrel{~}{m_\rho }^6+m_\eta ^4m_\varphi ^4+m_\varphi ^4\stackrel{~}{m_\rho }^4)\times `$ (18)
$`\left|{\displaystyle \frac{g_{\varphi \rho \pi }g_{\rho \eta \gamma }}{D_\rho (\stackrel{~}{m_\rho })}}\right|^2,`$ (19)
and

$`\stackrel{~}{m_\rho }^2=m_\eta ^2+{\displaystyle \frac{(m^2+m_\eta ^2m_\pi ^2)(m_\varphi ^2m^2)}{2m^2}}{\displaystyle \frac{(m_\varphi ^2m^2)x}{m}}p_{\pi \eta }`$ (20)
$`p_{\pi \eta }={\displaystyle \frac{\sqrt{(m^2(m_\eta m_\pi )^2)(m^2(m_\eta +m_\pi )^2)}}{2m}}.`$ (21)
Note that there is a misprint in Eq.(6) of Ref., which describes $`A_{back}(m,x)`$: the 7th term in the brackets ”$`+2m_\eta ^4\stackrel{~}{m_\rho }^4`$” should be replaced by ”$`+m_\eta ^4\stackrel{~}{m_\rho }^4`$”, as above in Eq.(19). Emphasize that all evaluations in Ref. were done with the correct formula.

The term of the interference between the signal and the background processes is written in the following way:

$$\frac{d\mathrm{\Gamma }_{int}(m)}{dm}=\frac{(m_\varphi ^2m^2)p_{\pi \eta }}{128\pi ^3m_\varphi ^3}_1^1𝑑xA_{int}(m,x),$$
(22)
where

$`A_{int}(m,x)={\displaystyle \frac{2}{3}}Re{\displaystyle }M_aM_B^{}={\displaystyle \frac{1}{3}}((m^2m_\varphi ^2)\stackrel{~}{m_\rho }^2+{\displaystyle \frac{m_\varphi ^2(\stackrel{~}{m_\rho }^2m_\eta ^2)^2}{m_\varphi ^2m^2}})\times `$ (23)
$`Re\{{\displaystyle \frac{e^{i\delta }g(m)g_{a_0K^+K^{}}g_{a_0\pi \eta }g_{\varphi \rho \pi }g_{\rho \eta \gamma }}{D_\rho ^{}(\stackrel{~}{m_\rho })D_{a_0}(m)}}\}.`$ (24)
Note that the phase $`\delta `$ isn’t taken into account in . The inverse propagator of the scalar meson R ($`a_0`$ in our case), is presented in Refs. :

$$D_R(m)=m_R^2m^2+\underset{ab}{}[Re\mathrm{\Pi }_R^{ab}(m_R^2)\mathrm{\Pi }_R^{ab}(m^2)],$$
(25)
where $`_{ab}[Re\mathrm{\Pi }_R^{ab}(m_R^2)\mathrm{\Pi }_R^{ab}(m^2)]=Re\mathrm{\Pi }_R(m_R^2)\mathrm{\Pi }_R(m^2)`$ takes into account the finite width corrections of the resonance which are the one loop contribution to the self-energy of the $`R`$ resonance from the two-particle intermediate $`ab`$ states.

For the pseudoscalar $`ab`$ mesons and $`m_am_b,mm_+`$ one has *** Note that in Ref. $`\mathrm{\Pi }_R^{ab}(m^2)`$ differs by a real constant from those determined in other enumerated works in the case of $`m_am_b`$, but obviously it has no effect on Eq.(25). :

$`\mathrm{\Pi }_R^{ab}(m^2)={\displaystyle \frac{g_{Rab}^2}{16\pi }}[{\displaystyle \frac{m_+m_{}}{\pi m^2}}\mathrm{ln}{\displaystyle \frac{m_b}{m_a}}+`$ (26)
$`+\rho _{ab}(i+{\displaystyle \frac{1}{\pi }}\mathrm{ln}{\displaystyle \frac{\sqrt{m^2m_{}^2}\sqrt{m^2m_+^2}}{\sqrt{m^2m_{}^2}+\sqrt{m^2m_+^2}}})]`$ (27)
For $`m_{}m<m_+`$

$`\mathrm{\Pi }_R^{ab}(m^2)={\displaystyle \frac{g_{Rab}^2}{16\pi }}[{\displaystyle \frac{m_+m_{}}{\pi m^2}}\mathrm{ln}{\displaystyle \frac{m_b}{m_a}}|\rho _{ab}(m)|+`$ (28)
$`+{\displaystyle \frac{2}{\pi }}|\rho _{ab}(m)|\mathrm{arctan}{\displaystyle \frac{\sqrt{m_+^2m^2}}{\sqrt{m^2m_{}^2}}}].`$ (29)
For $`m<m_{}`$

$`\mathrm{\Pi }_R^{ab}(m^2)={\displaystyle \frac{g_{Rab}^2}{16\pi }}[{\displaystyle \frac{m_+m_{}}{\pi m^2}}\mathrm{ln}{\displaystyle \frac{m_b}{m_a}}`$ (30)
$`{\displaystyle \frac{1}{\pi }}\rho _{ab}(m)\mathrm{ln}{\displaystyle \frac{\sqrt{m_+^2m^2}\sqrt{m_{}^2m^2}}{\sqrt{m_+^2m^2}+\sqrt{m_{}^2m^2}}}].`$ (31)
and

$$\rho _{ab}(m)=\sqrt{(1\frac{m_+^2}{m^2})(1\frac{m_{}^2}{m^2})},m_\pm =m_a\pm m_b$$
(32)
The constants $`g_{Rab}`$ are related to the width

$$\mathrm{\Gamma }(Rab,m)=\frac{g_{Rab}^2}{16\pi m}\rho _{ab}(m).$$
(33)
In our case we take into account intermediate states $`ab=\eta \pi ^0,K\overline{K}`$ and $`\eta ^{}\pi ^0`$:

$$\mathrm{\Pi }_{a_0}=\mathrm{\Pi }_{a_0}^{\eta \pi ^0}+\mathrm{\Pi }_{a_0}^{K^+K^{}}+\mathrm{\Pi }_{a_0}^{K^0\overline{K^0}}+\mathrm{\Pi }_{a_0}^{\eta ^{}\pi ^0},$$
(34)
$`g_{a_0K^+K^{}}=g_{a_0K^0\overline{K^0}}`$. Note that the $`\eta ^{}\pi ^0`$ contribution is of small importance due to the high threshold. Even fitting with $`|g_{a_0\eta ^{}\pi ^0}|=0`$ changes the results less than 10% of their errors. We set $`|g_{a_0\eta ^{}\pi ^0}|=|1.13g_{a_0K^+K^{}}|`$ according to the four-quark model, see , but this is practically the same as the 2-quark model prediction $`|g_{a_0\eta ^{}\pi ^0}|=|1.2g_{a_0K^+K^{}}|`$, see .

The inverse propagator of the $`\rho `$ meson has the following expression

$$D_\rho (m)=m_\rho ^2m^2im^2\frac{g_{\rho \pi \pi }^2}{48\pi }\left(1\frac{4m_\pi ^2}{m^2}\right)^{3/2}.$$
(35)
The coupling constants $`g_{\varphi K^+K^{}}=4.376\pm 0.074`$ and $`g_{\varphi \rho \pi }=0.814\pm 0.018`$ GeV-1 are taken from the new most precise measurement Ref.. Note that in Ref.  the value $`g_{\varphi K^+K^{}}=4.59`$ was obtained using the  data. The coupling constant $`g_{\rho \eta \gamma }=0.56\pm 0.05`$ GeV-1 is obtained from the data of Ref. with the help of the expression

$$\mathrm{\Gamma }(\rho \eta \gamma )=\frac{g_{\rho \eta \gamma }^2}{96\pi m_\rho ^3}(m_\rho ^2m_\eta ^2)^3.$$
(36)
## III Results

The KLOE data on the $`\varphi \eta \pi ^0\gamma `$ decay may be found in Table 5 of Ref. (see also Fig.1). Note that as in Refs., we don’t fit 1st,10th and 27th points of this table (cross points in the Fig. 1). Emphasize that the 10th (1.014 GeV) and 27th (1.019 GeV) points are obvious artifacts because the mass spectrum behaviour on the right slope of the resonance has the form the (photon energy)3 according to gauge invariance.

In the experiment the whole mass region ($`m_\eta +m_{\pi ^0},m_\varphi `$) is divided into some number of bins. Experimenters measure the average value $`\overline{B}_i`$ (”i” is the number of bin) of $`dBr(\varphi \eta \pi ^0\gamma )/dm`$ around each i-th bin:

$$\overline{B}_i=\frac{1}{m_{i+1}m_i}_{m_i}^{m_{i+1}}𝑑Br(\varphi \eta \pi ^0\gamma )/𝑑m,$$
(37)
In this case one should define $`\chi ^2`$ function as:

$$\chi ^2=\underset{i}{}\frac{(\overline{B}_i^{th}\overline{B}_i^{exp})^2}{\sigma _i^2},$$
(38)
where $`\overline{B}_i^{exp}`$ are the experimental results, $`\sigma _i`$ are the experimental errors, and

$$\overline{B}_i^{th}=\frac{1}{m_{i+1}m_i}_{m_i}^{m_{i+1}}𝑑Br^{th}(\varphi \eta \pi ^0\gamma )/𝑑m$$

($`dBr^{th}(\varphi \eta \pi ^0\gamma )/dm`$ is the theoretical curve).

The free parameters of the fit are $`m_{a_0},g_{a_0K^+K^{}}^2/4\pi `$, the phase $`\delta `$ (we assume it is constant) and the ratio $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}`$. The results are shown in Table 1 (fit 1) Note that fitting without averaging the theoretical curve $`(`$changing $`\overline{B}_i^{th}dBr^{th}(\varphi \eta \pi ^0\gamma )/dm|_{m=(m_{i+1}+m_i)/2})`$ results to worse $`\chi ^2/n.d.f.=28.8/20`$. The results in this case are consistent within errors with those obtained with averaging the theoretical curve..

The quality of the fit is good. The phase $`\delta `$ is consistent with zero, so we make a fit with $`\delta =0`$ (fit 2 in Table 1).

To check the correctness of treating the phase $`\delta `$ as a constant, we have done a fit with $`\delta `$ taken in the form $`\delta (m)=bp_{\eta \pi }(m)`$ (the phase of the elastic background in $`\eta \pi ^0`$ scattering may have such behaviour), and found out that the constant b=$`(2.8\pm 3.2)`$ GeV-1 is also consistent with zero. Change of the other values is not principal.

| Table 1. | | | |  |  |
| --- | --- | --- | --- | --- | --- |
| fit | $`m_{a_0}`$, MeV | $`\frac{g_{a_0K^+K^{}}^2}{4\pi }`$, GeV2 | $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}`$ | $`\delta `$, | $`\chi ^2/n.d.f.`$ |
| 1 | $`1003_{13}^{+32}`$ | $`0.82_{0.27}^{+0.81}`$ | $`1.06_{0.27}^{+0.20}`$ | $`27\pm 29`$ | 24.2/20 |
| 2 | $`995_8^{+22}`$ | $`0.65_{0.18}^{+0.42}`$ | $`1.17_{0.24}^{+0.17}`$ | 0 | 25.2/21 |
| 3 | $`994_8^{+22}`$ | $`0.62_{0.17}^{+0.4}`$ | $`1.21_{0.24}^{+0.17}`$ | $`21\pm 30`$ | 16.3/19 |
| 4 | $`992_7^{+14}`$ | $`0.55_{0.13}^{+0.27}`$ | $`1.26_{0.2}^{+0.16}`$ | $`0`$ | 16.9/20 |

Since the discrepancy between fits and the experimental point number 26 (0.999 GeV) in the Table 5 of Ref. (the box point in Fig.1) is about 3 standard deviations (i.e. this point may be an artifact also), we make another fit without this point (fit 3). The phase $`\delta `$ is again consistent with zero, so we make a fit without it (fit 4).

In Table 2 we present the results on the total branching ratio $`\text{Br}(\varphi (a_0\gamma +\rho \pi ^0)\eta \pi ^0\gamma )`$, the signal contribution $`\text{Br}(\varphi a_0\gamma \eta \pi ^0\gamma )`$, $`\mathrm{\Gamma }_{a_0\eta \pi ^0}\mathrm{\Gamma }(a_0\eta \pi ^0,m_{a_0})=g_{a_0\eta \pi }^2\rho _{\eta \pi ^0}(m_{a_0})/(16\pi m_{a_0})`$ and the ratio $`R=g_{f_0K^+K^{}}^2/g_{a_0K^+K^{}}^2`$. The last is obtained using the Ref. value $`g_{f_0K^+K^{}}^2/(4\pi )=2.79\pm 0.12`$ GeV2. The branching ratio of the background $`\text{Br}(\varphi \rho \pi ^0\eta \pi ^0\gamma )`$ accounts for $`(0.5\pm 0.1)\times 10^5`$.

| Table 2. | | | |  |
| --- | --- | --- | --- | --- |
| fit | $`\mathrm{\Gamma }_{a_0\eta \pi ^0}`$, MeV | R | $`10^5\times \text{Br}(\varphi (a_0\gamma +\rho \pi ^0)\eta \pi ^0\gamma )`$ | $`10^5\times \text{Br}(\varphi a_0\gamma \eta \pi ^0\gamma )`$ |
| 1 | $`153_{17}^{+22}`$ | $`3.4\pm 1.7`$ | $`7.6\pm 0.4`$ | $`7.3\pm 0.4`$ |
| 2 | $`148_{15}^{+17}`$ | $`4.3\pm 1.7`$ | $`7.6\pm 0.4`$ | $`7.1\pm 0.4`$ |
| 3 | $`149_{16}^{+19}`$ | $`4.5\pm 1.7`$ | $`7.6\pm 0.4`$ | $`7.2\pm 0.4`$ |
| 4 | $`146_{15}^{+17}`$ | $`5.0\pm 1.6`$ | $`7.6\pm 0.4`$ | $`7.1\pm 0.4`$ |

## IV Conclusion

Note that the obtained value of the ratio $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}`$ doesn’t contradict the first predictions based on the four-quark model of the $`a_0`$: $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}0.85`$ Note that the prediction $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}0.93`$, made in , was corrected in  . But even if $`g_{a_0\eta \pi }/g_{a_0K^+K^{}}`$ deviates from $`0.85`$, there is no tragedy, because those variant of the four-quark model is rather rough, it is considered as a guide.

For all fits the obtained value R differs from the  value $`R=7.0\pm 0.7`$. So the conclusion that the constant $`g_{a_0K^+K^{}}^2/(4\pi )`$ is small, obtained in  ($`g_{a_0K^+K^{}}^2/(4\pi )=0.4\pm 0.04`$ GeV2), is the result of the parameters restrictions, especially fixing $`m_{a_0}`$ at the PDG-2000 value 984.8 MeV. Note that a high $`a_0`$ mass is also needed to describe $`\gamma \gamma \eta \pi ^0`$ experiment, see .

There should be no confusion due to the large $`a_0`$ width. In the peripheral production of the $`a_0`$ (for example, in the reaction $`\pi ^{}p\eta \pi ^0n`$) the mass spectrum is given by the relation

$$\frac{dN_{\eta \pi ^0}}{dm}S_{per}(m)=\frac{2m^2}{\pi }\frac{\mathrm{\Gamma }(a_0\eta \pi ^0,m)}{|D_{a_0}(m)|^2}.$$
(39)
The effective (visible) width of this distribution is much less then the nominal width $`\mathrm{\Gamma }_{a_0\eta \pi ^0}`$. For example, for the fit 1 results (Table 1) the effective width is $`50`$ MeV, see Fig. 2.

As it is noted in Ref., there is no tragedy with the relation between branching ratios of $`a_0`$ and $`f_0`$ production in $`\varphi `$ radiative decays. The early predictions  are based on the one-loop mechanism $`\varphi K^+K^{}a_0\gamma \eta \pi ^0\gamma `$ and $`\varphi K^+K^{}f_0\gamma \pi \pi \gamma `$ at $`m_{a_0}=980`$ MeV, $`m_{f_0}=975`$ MeV and $`g_{a_0K^+K^{}}=g_{f_0K^+K^{}}`$§§§Emphasize that the isotopic invariance doesn’t require $`g_{a_0K^+K^{}}=g_{f_0K^+K^{}}`$., that leads to $`\text{Br}(\varphi a_0\gamma \eta \pi ^0\gamma )\text{Br}(\varphi f_0\gamma \pi \pi \gamma )`$. But it is shown in Ref. , that the relation between branching ratios of $`a_0`$ and $`f_0`$ production in $`\varphi `$ radiative decays essentially depends on a $`a_0f_0`$ mass splitting. This strong mass dependence is the result of gauge invariance, the (photon energy)3 law on the right slope of the resonance. Our present analysis confirms this conclusion. Note that a noticeable deviation from the naive four-quark model equality $`g_{a_0K^+K^{}}=g_{f_0K^+K^{}}`$ is not crucial. What is more important is the mechanism of the production of the $`a_0`$ and $`f_0`$ through the charged kaon loop, i.e. the four-quark transition. As is shown in Ref., this gives strong evidence in favor of the four-quark model of the $`a_0(f_0)`$.

Note that the constant $`g_{f_0K^+K^{}}^2/(4\pi )`$ also can differ a lot from those obtained in . The point is that the extraction of this constant is very model dependent. For example, fitting with taking into account the mixing of the resonances can decrease the value of $`g_{f_0K^+K^{}}^2/(4\pi )`$ considerably. For instance, fitting the data of  without mixing one has $`g_{f_0K^+K^{}}^2/(4\pi )=2.47_{0.51}^{+0.37}`$ GeV2 , while fitting with taking the mixing into account gives $`g_{f_0K^+K^{}}^2/(4\pi )=1.29\pm 0.017`$ GeV2 . Note also that in  the phase $`\delta _B`$ of the background is taken from the work , where it is obtained by the simultaneous fitting the $`m_{\pi ^0\pi ^0}`$ spectrum and the phase $`\delta _0`$ of the $`\pi \pi `$ scattering, taking into account the mixing of the resonances. In  the mixing isn’t taken into account, so additional phase dealing with it is omitted.

## V Acknowledgements

We thank C. Bini very much for providing the useful information, discussions and kind contacts. This work was supported in part by RFBR, Grant No 02-02-16061.


text:
# Curvaton Reheating in Non-oscillatory Inflationary Models

## Abstract

In non-oscillatory (NO) inflationary models, the reheating mechanism was usually based on gravitational particle production or the mechanism of instant preheating. In this paper we introduce the curvaton mechanism into NO models to reheat the universe and generate the curvature perturbations. Specifically we consider the Peebles-Vilenkin quintessential inflation model, where the reheating temperature can be extended from $`1`$MeV to $`10^{13}`$GeV.

PACS number(s): 98.80.Cq

Inflarionary universe guth  has solved many problems of the Standard Hot Big-Bang theory, e.g., flatness problem, horizon problem, etc. In addition, it has provided a causal interpretation of the origin of the observed anisotropy of the cosmic microwave background (CMB) and the distribution of large scale structure (LSS). In the usual inflation models, the acceleration is driven by the potential of a scalar field $`\varphi `$ (inflaton) and its quantum fluctuations in this epoch generate the density perturbations seeding the structure formations at late time. To date, the accumulating observational data, especially those from the CMB observation of WMAP satellite wmap1 ; wmap2 ; wmap3  indicate the power spectrum of the primordial density perturbations is nearly scale-invariant, adiabatic, and Gaussian–just as predicted by the single-field inflation in the context of “slow-roll”. Generally, after inflation the inflaton field will oscillate about the minimum of its potential and eventually decay to produce almost all elementary particles populating the universe. This process is called reheating.

However, there exist some inflationary models no ; ford ; vilenkin ; joyce , named NO models by the authors of Ref. kofman , in which the minimum of the inflaton potential exist in the infinity. In these models the inflaton field would not oscillate after inflation and the standard reheating mechanism could not work. Conventionally one turns to gravitational particle production ford . However, the authors of Ref. kofman  pointed out that this mechanism is very inefficient. For the resolution to this problem, they proposed to use “instant preheating” felder  by introducing an interaction $`g^2\varphi ^2\chi ^2`$ of the inflaon field $`\varphi `$ with another scalar field $`\chi `$.

Furthermore, in these models since the inflaton field can not decay away efficiently, its energy density must decrease extremely quickly ($`\mathrm{e}.\mathrm{g}.,\rho _\varphi a^6`$) soon after the end of inflation so that it cannot dominate the universe again too early. Otherwise, the successful achievements of big-bang nucleosynthesis (BBN), CMB and structure formations would be spoiled. This rapid diminishing phase was called “kination” in Ref. joyce  or “deflation” in Ref. no . It requires a steep region in the potential connecting to the flat plateau for inflation. This often leads to the potential being strongly curved near the end of inflation and tend to result in a power spectrum too far from scale-invariance in some NO models. This problem is called $`\eta `$-problem in the literature dimo ; liddle2 , because in these models the predicted spectral indices

$$n_s=1+2\eta 6ϵ,$$
(1)
are in conflict with current observational value $`n_s=0.93\pm 0.03`$ wmap1 ; wmap2  mainly due to a large $`|\eta |`$. In Eq. (1) the slow-roll parameters are defined as

$$ϵ\frac{M_{pl}^2}{16\pi }(\frac{V^{}}{V})^2,\eta \frac{M_{pl}^2}{8\pi }\frac{V^{\prime \prime }}{V},$$
(2)
where $`M_{pl}=1.22\times 10^{19}`$GeV is the Planck mass and prime denotes the derivative with respect to $`\varphi `$. In addition, in a concrete NO model, the so-called original “quintessentail inflation” vilenkin  (which will be quoted as the example in present paper), the inflation epoch is the same as chaotic inflation with potential $`V=\lambda \varphi ^4`$. It predicts a large gravitational wave amplitude, disfavored by the WMAP data wmap3 .

In this paper we put forward another reheating mechanism for NO models by introducing the curvaton mechanism lyth  which received many attentions recently in the literature dimo ; liddle2 ; enqv ; liddle ; wands ; moroi ; enqvist , and the problems mentioned above can be ameliorated. In our scenario, the reheating mechanism is based on the oscillations and decays of another scalar field $`\sigma `$ (the curvaton field) which is sub-dominant during inflation and has no interactions with inflaton except the gravitational coupling. This process is as same as that of the standard reheating in the usual inflation models. We will show by an example that the reheating temperature in our case can be as high as $`10^{13}`$GeV. In addition, in the curvaton mechanism the curvature (adiabatic) perturbations are not originated from the fluctuations of the inflton field, but instead from those of curvaton lyth . This happens in two steps. First, the quantum fluctuations in the curvaton field during inflation are converted into classical perturbations, corresponding to isocurvature perturbations when cosmological scales leave the horizon. Then, after inflation when the energy density of the curvaton field becomes significant, the isocurvature perturbations in it are converted into sizable curvature ones mollerach . The curvature perturbations become pure when the curvaton field or the radiation it decays into begins to dominate the universe. Because the curvaton is very light during inflation, the final spectrum of the curvature perturbations is nearly flat and naturally consistent with current observations. At the same time, in curvaton scenario the fluctuations in inflaton field is negligible, so the energy scale of inflation is relatively lower and the gravitational wave amplitude is very small. Existing observations have detected no gravitational waves, hence those inflation models which usually predicted unlikely large tensor perturbations based on the standard paradigm of inflaton-generating curvature perturbations, like $`\lambda \varphi ^4`$ (and the example NO model cited below) wmap3 , would become viable in the curvaton scenario. The primordial density perturbations in curvaton scenario have been studied specifically in Ref. wands  and its effects on CMB analyzed by the authors of Ref. moroi ; lewis . Also, the curvaton mechanism has been used to liberate some inflation models liddle ; dimo . In addition, there were some motivations for the curvaton model from particle physics lyth ; enqvist .

For a specific presentation, we consider the quintessential inflation model given by Peebles and Vilenkin vilenkin . In this model, the inflaton field will survive to now to drive the present accelerated expansion suggested by recent measurements of type Ia supernova super . So, it can not decay to other particles. The potential of the infaton field is:

$`V(\varphi )`$ $`=`$ $`\lambda (\varphi ^4+M^4)\mathrm{for}\varphi <0`$ (3)
$`=`$ $`{\displaystyle \frac{\lambda M^8}{\varphi ^4+M^4}}\mathrm{for}\varphi 0.`$
The energy scale $`M`$ is assumed very small compared with the Planck mass, the requirement of $`\varphi `$ driving the present accelerated expansion gives that $`\varphi _0>M_{pl}`$ and $`M>10^5`$GeV, the subscript $`0`$ represents present value. The inflation which happened at $`\varphi M`$ is chaotic type with $`V=\lambda \varphi ^4`$. The slow-roll parameters are

$$ϵ=\frac{M_{pl}^2}{\pi \varphi ^2},\eta =\frac{3M_{pl}^2}{2\pi \varphi ^2}.$$
(4)
Inflation ends when $`ϵ=1`$, this corresponds to $`\varphi _{end}0.56M_{pl}`$. The e-folding number between the horizon-exit of the comoving mode and the termination of inflation can be estimated by

$$N_{t_N}^{t_{end}}H𝑑t\frac{\pi \varphi _N^2}{M_{pl}^2}1.$$
(5)
Hence,

$$\varphi _N\sqrt{\frac{N+1}{\pi }}M_{pl},$$
(6)
and

$$ϵ(\varphi _N)\frac{1}{N+1}.$$
(7)
The power spectra of generated curvature and tensor perturbations by inflaton are respectively book ,

$`𝒫_\zeta ={\displaystyle \frac{8V}{3M_{pl}^4ϵ}},`$
$`𝒫_h={\displaystyle \frac{128V}{3M_{pl}^4}}.`$ (8)
If the cosmic density perturbations of order $`10^{10}`$ rely completely on this curvature perturbation, the tensor to curvature ratio $`r𝒫_h/𝒫_\zeta =16ϵ`$ will be too large (for $`N=50`$, $`r=0.32`$) to be compatible with WMAP data at $`3\sigma `$ wmap3 .

This difficulty can be alleviated in the curvaton mechanism. In this mechanism the curvature perturbations produced by the inflaton is supposed to be negligible in contrast with COBE data $`𝒫^{1/2}4.8\times 10^5`$. This leads to

$`\lambda <10^{14},`$
$`H_{in}<10^6M_{pl},`$ (9)
where the subscript $`in`$ denotes the value taken by the corresponding variable in the inflationary epoch. So, the gravitational wave amplitude, proportional to $`H_{in}/M_{pl}`$, is negligibly small. In the original paper of curvaton scenario, only a part of the components of the universe come from the decays of the curvaton field. In our case, since inflaton could not decay, we can assume that the whole population of the universe is produced by curvaton. We will show that this process can be very efficient.

We choose the quadratic potential for the curvaton field, $`V(\sigma )=m^2\sigma ^2/2`$, it has a minimum in the zero point of $`\sigma `$ and $`mH_{in}`$.

We numerically studied the evolution of the energy densities in the inflton field and in the curvaton field. As demonstrated in Fig. 1, in the era of inflation, the curvaton field slowly rolled towards the minimum along its potential and its energy density changed little. Inflation ended at $`a=a_{end}`$ when $`\varphi 0.5M_{pl}`$, $`a`$ is the scale factor. Then the universe entered into the kination period dominated by the kinetic energy of inflaton, during which:

$`\rho _\varphi a^6,Ha^3.`$ (10)
The curvaton field held slow-rolling at the beginning of kination because of $`mH`$, and its energy density $`\rho _\sigma `$ decreased very slowly. Then it began to oscillate about its minimum at $`a=a_{osc}`$ when $`H`$ decreased below $`m`$. To avoid another inflation stage driven by the curvaton field, the universe must still be dominated by $`\rho _\varphi `$ at this point, this requires

$$\rho _{\sigma osc}<\rho _{\varphi osc}=\frac{3H_{osc}^2M_{pl}^2}{8\pi }.$$
(11)
From

$$H_{osc}m,\rho _{\sigma osc}m^2\sigma _{osc}^2,$$
(12)
we can estimate that:

$$\sigma _{osc}^2<\frac{3M_{pl}^2}{8\pi }.$$
(13)
After the curvaton field began oscillating, $`\rho _\varphi `$ continued to decrease as $`a^6`$ towards the value of present energy density of dark energy and subsequently froze almost like a cosmological constant when $`aa_f`$. $`\rho _\sigma `$ red-shifted as $`a^3`$ like that of matter when averaged over many oscillations until it decayed into radiation at $`a=a_{rh}`$, the subscript $`rh`$ means reheating. This happened when $`H=H_{rh}\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }`$ is the decay rate of the curvaton field (we have assumed instant reheating). It can occur after the moment when the curvaton field dominated the universe ($`a_{rh}a_{eq}`$) or at the epoch when the curvaton field was still sub-dominant ($`a_{rh}<a_{eq}`$).

In the first case, as shown in Fig. 1, this requires:

$$\rho _{\varphi rh}\rho _{\sigma rh}.$$
(14)
Since $`\rho _\varphi a^6`$, $`\rho _\sigma a^3`$ and

$$H_{osc}/H_{eq}=(a_{eq}/a_{osc})^3,H_{eq}^2/H_{rh}^2=(a_{rh}/a_{eq})^3,$$
(15)
we get

$$\frac{\mathrm{\Gamma }^2}{mH_{eq}}\frac{8\pi \sigma _{osc}^2}{3M_{pl}^2}.$$
(16)
In above we have used Eq. (12). From $`\rho _{\sigma eq}=\rho _{\varphi eq}`$ we have

$$H_{eq}=\frac{8\pi m\sigma _{osc}^2}{3M_{pl}^2},$$
(17)
so the constraint is (with Eq. (13)):

$$\frac{\mathrm{\Gamma }}{m}\frac{8\pi \sigma _{osc}^2}{3M_{pl}^2}<1.$$
(18)
The reheating temperature can be estimated from:

$$\rho _{\sigma rh}=\frac{3M_{pl}^2}{8\pi }\mathrm{\Gamma }^2\rho _{rad}=\frac{\pi ^2}{30}g_{}T_{rh}^4,$$
(19)
where $`g_{}`$ counts the total number of degrees of freedom of the relativistic particles, $`g_{}=10.75`$ at temperature $`T1`$MeV, and for $`T100`$GeV, $`g_{}𝒪(100)`$. The reheating process must be completed before BBN, so in this case, the reheating temperature is approximately

$$T_{rh}0.78g_{}^{1/4}\sqrt{M_{pl}\mathrm{\Gamma }}\sqrt{m\sigma _{osc}^2/M_{pl}}.$$
(20)
Similarly, in the latter case, when the curvaton field decayed when it was still sub-dominant but later than oscillation ($`a_{osc}a_{rh}<a_{eq}`$), we have (for general, we assume $`\sigma _{osc}>0`$)

$`{\displaystyle \frac{8\pi \sigma _{osc}^2}{3M_{pl}^2}}<{\displaystyle \frac{\mathrm{\Gamma }}{m}}<1,`$
$`\sqrt{m\sigma _{osc}^2/M_{pl}}<T_{rh}<0.4\sqrt{m\sigma _{osc}}.`$ (21)
In this case, the radiation dominated era began with $`T_{eq}(m\sigma ^2)^{3/4}/(\mathrm{\Gamma }^{1/4}M_{pl})`$ which is in the range

$$\sqrt{m\sigma _{osc}^3}/M_{pl}<T_{eq}<0.6\sqrt{m\sigma _{osc}^2/M_{pl}}.$$
(22)
For evaluation of $`T_{rh}`$, one has to consider the constraints on the mass of $`\sigma `$ field and its initial value when it began oscillating, $`\sigma _{osc}`$. Based on the curvaton mechanism lyth , the resulted power spectrum of the final curvature perturbations and its index are given by

$`𝒫_\zeta ^{1/2}={\displaystyle \frac{H_{in}}{3\pi \sigma _{in}}}0.1{\displaystyle \frac{H_{in}}{\sigma _{in}}},`$ (23)
$`n_s1+d\mathrm{ln}𝒫_\zeta /d\mathrm{ln}k=12ϵ+2m^2/(3H_{in}^2).`$ (24)
One can see that the spectrum is nearly flat. The curvaton field rolled very slowly before oscillation, we can evaluate that $`\sigma _{osc}\sigma _{in}`$. Combined these results with the COBE normalization and the limit $`H_{in}<10^6M_{pl}`$, one has

$$\sigma _{osc}<10^3M_{pl}10^{16}\mathrm{GeV}.$$
(25)
Substituting Eq. (25) into Eq. (20) and (Curvaton Reheating in Non-oscillatory Inflationary Models), and considering $`mH_{in}`$ (we choose $`m10^7M_{pl}10^{12}`$GeV for an evaluation), in both cases we have approximately

$$T_{rh}10^{13}\mathrm{GeV}.$$
(26)
To ensure above analytical estimations we have made a simple numerical calculation. We set $`\lambda =10^{15}`$, $`M=4.7\times 10^{14}M_{pl}`$, $`m_\sigma =10^8M_{pl}`$ and $`\sigma _{in}=3.27\times 10^3M_{pl}`$. For curvaton decays at $`a=a_{eq}`$, one gets $`T_{rh}1.8\times 10^{12}`$GeV. Assuming that curvaton decays into two fermions $`\sigma f\overline{f}`$, the decay width takes the form $`\mathrm{\Gamma }=g^2m_\sigma /8\pi `$. We find for $`\sigma `$ decays at $`a=a_{eq}`$ one requires $`g0.05`$, which is quite resonable. When the curvaton decays at $`a<a_{eq}`$, $`g<1`$ gives $`T_{rh}<8\times 10^{12}`$GeV. The BBN constraint that $`T_{rh}>1`$MeV can also be easily satisfied, as from Fig. 1.

In conclusion, we have introduced the curvaton mechanism into the non-oscillatory inflationary models as a possible reheating mechanism and did specific studies in the quintessential inflation model proposed by Peebles and Vilenkin. In our scenario, the reheating process is not attributed to the gravitational particle productions, but due to the late-decaying of the curvaton field. This mechanism, as same as the standard reheating mechanism which based on the oscillations and decays of the scalar field (inflaton) in the usual chaotic inflation models, can be efficient. Our calculations showed that the reheating temperature can be as high as $`10^{13}`$GeV, it is high enough for the quintessential baryogenesis li ; trodden . At the same time, in addition to reheat the universe, curvaton has been endowed the role of generating curvature perturbations as in the original paper lyth . By introducing curvaton mechanism into the NO models, the predicted power spectra of curvature perturbations are nearly scale-invariance and the tensor perturbations are negligibly small, consistent with current observations. These properties can ameliorate other defects owned by some NO models. For example, in Peebles-Vilenkin model we quoted in this paper, the inflationary period is as same as the chaotic inflation $`V=\lambda \varphi ^4`$, it predicts unlikely large gravitational waves according to the usual paradigm of inflaton-generating curvature perturbations wmap3 . For some other NO models, the density perturbations are far from scale-invariance dimo ; liddle2 , in conflict with observations if the cosmic inhomogeneities are only due to the fluctuations of the inflaton field. Our scenario provides a picture that the particle productions and the cosmic density perturbations are both due to the curvaton field, this made it different from the original curvaton scenario (in which only a part of the components of the universe come from curvaton) and other reheating mechanisms that had been brought forward for NO inflationary models. Our work has also provided a new way to the reheating problem of tachyon inflation tachyon .

While this work was in progress a related paper dimo  appeared. In Ref. dimo  the author mainly dealt with qunintessential inflation model building and parameter restriction and our paper is focusing on the curvaton reheating mechanism.

Acknowledgments: We would like to thank Prof. Xinmin Zhang and Dr. Yun-Song Piao for useful discussions. We also thank Prof. David Lyth and Prof. David Wands for correspondences. This work is supported in part by National Natural Science Foundation of China and by Ministry of Science and Technology of China under Grant No. NKBRSF G19990754.


text:
# Teleportation of a Bose-Einstein condensate state by controlled elastic collisions

## I Introduction

Teleportation of quantum states, proposed by Bennett et al. , was first realised for light polarization states , owing to the possibility of generating nonlocal entanglement between parties (the quantum channel) in this system. Although many proposals and experimental realizations of nonlocal entanglement of massive particles (atoms and ions) exist , up to now there has been no experimental evidence of the teleportation of massive particles state . Indeed, non-linear interactions (a valuable resource for deterministic generation of entanglement) are always present in many-particle systems. An interesting question then arises - to what extent can the teleportation protocol be applied to massive many-particle systems? A strong candidate for massive particle state teleportation is the condensate state of a matter field, where non-linear interactions appear quite naturally as elastic collisions .

In this paper we propose an experimental protocol for teleportation  of mode states of an atomic Bose-Einstein condensate (BEC) trapped in an optical lattice potential, by using controllable elastic collisions and Josephson coupling  between modes for both the quantum channel formation and measurement stage. Elastic collisions are a fundamental resource for both the formation of the entanglement and the parity operations needed to correct the teleported state. For the measurement stage, we propose a balanced homodyne detection of the BEC modes quadratures, valid for a small condensate. By measuring the difference of population in two condensate modes (the central mode and a reference mode) interacting via Josephson coupling the central mode quadrature is determined.

This paper is organized as follows. In Sec. II we present the model for three interacting condensate modes trapped in a three-well potential. In Sec. III we propose the teleportation protocol using controlled collisions as a means to generate entanglement. In Sec. IV we turn to the measurement stage, presenting a scheme of homodyne detection of the BEC phase. In Sec. V we present the operations of parity and displacement needed in order to correct to the required state, the state of the condensate mode at the receiver station. In Sec. VI we present a physical implementation of the controlled collision on optical lattices. In Sec. VII a discussion encloses the paper.

## II model

The observation of BECs of diluted trapped neutral atoms  and the recent achievement of all-optical confinement of a BEC  and condensation on a microelectronic chip  have stimulated a large research program on BEC of diluted neutral atom gases. Of particular interest is the study of a BEC in a confining potential. In  the dynamics of a BEC in a double-well potential was modelled. In that model, coherent oscillations due to tunnelling (Josephson-like coupling)  between the two wells are suppressed when the number of atoms exceeds a critical value (self-trapping). In fact, the barrier separation between the two wells has a central role as it determines whether the Josephson coupling between modes is negligible in contrast to cross-collisions, when the wave-functions of the two modes considerably overlap. A dynamic process can be envisaged where the two-well barrier is lowered and raised back adiabatically, such that the elastic collisions leave the two modes in an entangled state  \- a hallmark of teleportation protocols. This is the mechanism we focus on here. We extensively use controlled collisions between condensate modes, by adiabatically turning off and on the potential barriers in an optical lattice potential, in order to teleport the state of one condensate mode, to another mode, located inside the trap, but delocalised from the first one. This sequential process is depicted in Fig. 1. Initially we describe the general approach for entanglement generation and measurement and at the end we give a specific but clarifying physical implementation of the time varying potential, for optical lattices.

Extending the model in , we consider a BEC trapped in a symmetric three-well single-particle potential $`V(r)`$ with minima at $`r_1`$, $`r_2`$ and $`r_3`$ disposed along the $`z`$ axis. We assume that the three lowest states of the potential are closely spaced and well separated from its higher levels, and that many-particle interactions do not significantly change this situation, to allow a three-mode approximation. The potential expanded around each minimum is

$$V(𝐫)=\stackrel{~}{V}^{(2)}(𝐫𝐫_j)+\mathrm{}j=1,2,3,$$
(1)
where $`\stackrel{~}{V}^{(2)}(𝐫𝐫_j)`$ is the parabolic approximation to the potential in the vicinity of each minimum. The normalized single-particle ground-state $`u_0(𝐫)`$ of the local potential $`\stackrel{~}{V}^2(𝐫)`$, with energy $`E_0`$, defines the local mode solutions of the individual wells. If the position uncertainty in the state $`u_0(𝐫)`$ is much less than the separation of the minima of the global potential, the overlap between the modes of each well, $`ϵ`$, is much less than unity and the modes are approximately orthogonal. The many-body Hamiltonian describing an atomic BEC in this potential is

$`H`$ $`=`$ $`{\displaystyle d^3x\psi ^{}(𝐫)\left(\frac{\mathrm{}}{2m}^2+V(𝐫)\right)\psi (𝐫)}`$ (3)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{4\pi a_s\mathrm{}^2}{m}}{\displaystyle d^3r\psi ^{}(𝐫)\psi ^{}(𝐫)\psi (𝐫)\psi (𝐫)},`$
where $`m`$ is the atomic mass, $`U_0=4\pi \mathrm{}^2a/m`$ measures the strength of the two-body interaction, $`a`$ being the s-wave scattering length, $`\psi (𝐫,t)`$ and $`\psi ^{}(𝐫,t)`$ are the Heisenberg picture field operators, which annihilate and create atoms at position $`𝐫`$, and normal ordering has been used. In the three-mode approximation the field operators are expanded in terms of the local modes and the Heisenberg picture annihilation and creation operators read as

$$c_j(t)=d^3𝐫u_j^{}(𝐫)\psi (𝐫,t)$$
(4)
so that $`[c_j,c_k]=\delta _{jk}`$ to order $`ϵ^0`$. With this prescription, and retaining terms up to order $`ϵ`$, one arrives at two distinct regimes: (i) If the potential wells are well separated only self-collision terms are important and the many-body Hamiltonian reduces to

$`H_1`$ $`=`$ $`E_0(c_1^{}c_1+c_2^{}c_2+c_3^{}c_3)`$ (6)
$`+\mathrm{}\kappa [(c_1^{})^2c_1^2+(c_2^{})^2c_2^2+(c_3^{})^2c_3^2],`$
where $`\kappa =U_0/2\mathrm{}V_{eff}`$, and $`V_{eff}^1=d^3𝐫|u_0(𝐫)|^4`$ is the effective mode volume of each well. In such a situation, no cross-collisions or Josephson tunnelling occur. (ii) If the potential wells are not well separated Josephson tunnelling  between neighbour wells occurs and in the diluted atomic gas regime it prevails over cross-collisions. The many-body Hamiltonian then becomes

$`H_2`$ $`=`$ $`E_0(c_1^{}c_1+c_2^{}c_2+c_3^{}c_3)+\mathrm{}{\displaystyle \frac{\mathrm{\Omega }}{2}}(c_1^{}c_2+c_2^{}c_3+H.c.)`$ (8)
$`+\mathrm{}\kappa [(c_1^{})^2c_1^2+(c_2^{})^2c_2^2+(c_3^{})^2c_3^2],`$
where $`\mathrm{\Omega }=2/\mathrm{}`$ is the tunnelling frequency between two minima, with

$$=d^3𝐫u_i^{}(𝐫)[V(𝐫)\stackrel{~}{V}^{(2)}(𝐫𝐫_i)]u_{i+1}(𝐫),i=1,2.$$
(9)
We assume that in equilibrium state, regime (i) rules out, and the modes can be treated independently. If the barrier separating modes 2 and 3 is lowered adiabatically ($`\frac{dV}{dt}(E^{}E_0)/\mathrm{}`$), in order to avoid transitions to other states (of energy $`E^{}`$), the two respective modes overlap and strong cross collision occurs as $`V_c=2\mathrm{}\kappa c_2^{}c_2c_3^{}c_3`$. When the barrier is raised back adiabatically, bringing the system to the equilibrium regime, the two modes split again, but now having a non-local entanglement generated by the cross collision term, over the time the two modes overlapped. In such a case a non-local quantum channel would be formed between condensate modes 2 and 3  as we describe in next section.

## III Teleportation Protocol

The following protocol is more efficient if the condensate modes are initially prepared in coherent states, but we may assume a general initial state expanding it in the coherent state basis, $`|\psi _a_2=d^2\alpha a_\alpha |\alpha _2`$ and $`|\psi _b_3=d^2\beta b_\beta |\beta _3`$, for modes 2 and 3 respectively. The dynamics governed by $`H_1`$, together with the cross-collision $`V_c`$, gives for modes 2 and 3 alone

$`|\psi (t)`$ $`=`$ $`{\displaystyle d^2\alpha d^2\beta a_\alpha b_\beta e^{(|\alpha |^2+|\beta |^2)/2}}`$ (12)
$`\times {\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\sqrt{m!n!}}}\left(\alpha e^{i(E_0/\mathrm{}\kappa )t}e^{i\kappa mt}\right)^m`$
$`\times \left(\beta e^{i(E_0/\mathrm{}\kappa )t}e^{i\kappa (m+n)t}\right)^n|m_2|n_3,`$
which, for $`t=\pi /2\kappa `$, turns out to be the entangled state given by

$`|\mathrm{\Phi }(\pi /\kappa )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^2\alpha d^2\beta a_\alpha b_\beta [(1i)|\alpha e^{i\varphi }_2|\beta e^{i\varphi }_3`$ (14)
$`+(1+i)|\alpha e^{i\varphi }_2|\beta e^{i\varphi }_3]`$
where $`\varphi =(E_0/\mathrm{}\kappa )/2\kappa `$. Choosing properly the frequency of the modes, $`(E_0/\mathrm{}\kappa )`$, a set of approximately orthonormal states $`|\mathrm{\Phi }_j`$ can be generated for $`E_0/\mathrm{}=(j+1)\kappa `$, $`j=0,1,2,3`$, respectively as

$`|\mathrm{\Phi }_j`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^2\alpha d^2\beta a_\alpha b_\beta [(1i)|(i)^j\alpha _2|(i)^j\beta _3`$ (16)
$`+(1+i)|(i)^j\alpha _2|(i)^j\beta _3].`$
From now on we suppose that the condition for the generation of $`|\mathrm{\Phi }_0`$, i.e., $`E_0/\mathrm{}=\kappa `$, is met.

A simple teleportation protocol  based on homodyne measurement of the condensate modes phases can be performed. Let us say the condensate mode 1, called hereafter the target mode, is prepared in an unknown superposition of the specific form

$$|\psi _T=d^2\gamma c_\gamma \left(A|\gamma +B|\gamma \right)$$
(17)
where $`A`$ and $`B`$ are constants respecting normalization conditions. Now if the modes 1 and 2 are made to collide the whole condensate state is left as

$`{\displaystyle \frac{1}{2}}{\displaystyle }d^2\alpha d^2\beta d^\gamma a_\alpha b_\beta c_\gamma \{i|\gamma |\alpha (A|\beta B|\beta )`$ (18)
$`+|\gamma |\alpha (A|\beta +B|\beta )+i|\gamma |\alpha (A|\beta B|\beta )`$ (19)
$`+|\gamma |\alpha (A|\beta +B|\beta )\},`$ (20)
a three-partite entangled state composed of four elements. We can distinguish each element by the phase of modes 2 and target. Notice that although the protocol is encoded in continuous variables states, the protocol itself is discrete as only four equiprobable outcomes are possible. Thus after the distinction of the target and mode 2 phases, only two bits of classical information have to be sent to the mode 3 at the receiver station.

## IV Homodyne detection of BEC phase

If a joint measurement on both modes 1 and 2 can be envisaged to distinguish the phase of each condensate mode, the complete Bell state measurement is realised. Tomographic reconstruction of the condensate state would allow the distinction between the two different phases. Such an approach is based on optical homodyne tomography , using an arrangement composed of an atomic beam splitter and an ideal atom counter. Here instead, we describe an alternative scheme for phase determination similar the optical balanced homodyne measurement, where the Josephson coupling  plays the role of an atomic beam splitter. Consider two condensate modes separated by a barrier, as in . A two-mode approximation is assumed. Neglecting cross-collision terms (once the overlap of condensate wave-functions is negligible) the Hamiltonian for the two modes is

$`H`$ $`=`$ $`E_0(c^{}c+b^{}b)+{\displaystyle \frac{\mathrm{}\mathrm{\Omega }}{2}}(c^{}b+b^{}c)`$ (22)
$`+\mathrm{}\kappa [(c^{})^2c^2+(b^{})^2b^2],`$
Defining new operators as $`S_x=\frac{1}{2N}(c^{}cb^{}b)`$, $`S_y=\frac{i}{2N}(c^{}bcb^{})`$, $`S_z=\frac{1}{2N}(c^{}b+cb^{})`$, where $`N=c^{}c+b^{}b`$, the equations for the evolution of these operators are

$`\dot{S}_x`$ $`=`$ $`\mathrm{\Omega }S_y,`$ (23)
$`\dot{S}_y`$ $`=`$ $`\mathrm{\Omega }S_x2iϵ\mathrm{\Omega }S_y4ϵ\mathrm{\Omega }NS_xS_z,`$ (24)
$`\dot{S}_z`$ $`=`$ $`2iϵ\mathrm{\Omega }S_z+4ϵ\mathrm{\Omega }NS_xS_y,`$ (25)
where $`ϵ=\frac{\kappa }{\mathrm{\Omega }}1`$. A semiclassical solution for $`S_x`$ is given up to first-order in $`ϵ`$ by

$`S_x(t)`$ $`=`$ $`[S_x(0)+ϵt(2Nz_0y_0ix_0)]\mathrm{cos}\mathrm{\Omega }t`$ (27)
$`[S_y(0)ϵt(2Nz_0x_0+iy_0)]\mathrm{sin}\mathrm{\Omega }t,`$
valid only for $`ϵN1`$, i.e., for $`\kappa \mathrm{\Omega }`$, a small number of particles and also for a short time. To derive this solution the above operators are expanded as $`S_x=_nϵ_nx_n`$, $`S_y=_nϵ_ny_n`$ and $`S_z=_nϵ_nz_n`$. Assuming initially an equal number of atoms in both wells, the solution (27) simplifies to

$`S_x(t)`$ $`=`$ $`S_y(0)\mathrm{sin}\mathrm{\Omega }t+2ϵtNz_0y_0\mathrm{cos}\mathrm{\Omega }t+iϵty_0\mathrm{sin}\mathrm{\Omega }t.`$ (28)
Consider the following semiclassical picture for the operator $`S_y`$

$`S_y={\displaystyle \frac{i}{2}}|\beta |(c^{}e^{i\theta }ce^{i\theta })=|\beta |X_{\theta \pi /2},`$ (29)
where the mode B was prepared in a coherent state given by $`\beta =|\beta |e^{i\theta }`$. It is easy to observe then that for $`ϵN1`$, at $`t=\pi /2\mathrm{\Omega }`$, the Eq.(28) gives the well-known result for balanced homodyne detection, plus a small correction

$$S_x(\pi /2\mathrm{\Omega })=|\beta |X_{\theta \pi /2}+i\frac{\pi ϵ}{2\mathrm{\Omega }}y_0,$$
(30)
i.e., the difference between the numbers of atoms in the two wells determines the quadrature phase of one of the matter fields. This method is, however, sensitive to the exact determination of the reference phase $`\theta `$, which for condensates is a central problem . Here it is simply assumed that the reference phase can be determined by the experimentalist. With such a quadrature matter field measurement at hand it is possible to distinguish between coherent states like $`|\alpha `$ and $`|\alpha `$, which is the necessary resource to apply to both target and mode 2 and thus distinguish between the many states of the superposition of Eq.(18). Remark that the requirement for a small number of particles ($`N\mathrm{\Omega }/\kappa `$) avoids the regime of self-trapping, as observed in , when the approximate solution, Eq. (27), is no longer valid. This requirement imposes a severe limitation of this detection method to the “size” of the BEC to be teleported.

## V Receiver operations

The two bits of classical information obtained as described above are transferred to the receiver mode 3. Now one has to apply the operations needed to transform the condensate state mode 3 in the receiver station into the required state. Depending on the results of the joint homodyne measurement described above, the condensate mode 3 is left in one of the following states,

$`{\displaystyle }d^2\beta b_\beta (A`$ $`|\beta +B|\beta `$ $`)`$ (31)
$`{\displaystyle }d^2\beta b_\beta (A`$ $`|\beta B|\beta `$ $`)`$ (32)
$`{\displaystyle }d^2\beta b_\beta (A`$ $`|\beta +B|\beta `$ $`)`$ (33)
$`{\displaystyle }d^2\beta b_\beta (A`$ $`|\beta B|\beta `$ $`)`$ (34)
In order to transform the condensate states (32), (33), and (34) into the required state (31), operations of parity and displacement , in principle, can be realised, as described below.

### A Parity

A parity operation involving only atomic systems may be envisaged for a two-species BEC with the cross-collision strenght given by $`\lambda =U_{AB}/2\mathrm{}V_{eff}`$, where $`\lambda \kappa =U_{AA}/2\mathrm{}V_{eff}=U_{BB}/2\mathrm{}V_{eff}`$ . Supposing the central condensate initially in the state (33) and we require to transform it to (31). The auxiliary condensate is prepared in an arbitrary state expanded also in the coherent basis $`|\psi _a=d^2\xi a_\xi |\xi `$. Due to cross-collision the state of this system at time $`t=2\pi /\kappa `$ is

$`|\mathrm{\Phi }_{ac}`$ $`=`$ $`{\displaystyle d^2\beta d^2\xi };b_\beta a_\xi {\displaystyle \underset{n}{}}{\displaystyle \frac{e^{\frac{1}{2}|\xi |^2}}{\sqrt{n!}}}(\xi e^{i2\pi \frac{(E_0/\mathrm{}\kappa )}{\kappa }})^n|n_a`$ (36)
$`(A|\beta e^{i2\pi \frac{(E_0/\mathrm{}\kappa +n\lambda )}{\kappa }}+B|\beta e^{i2\pi \frac{(E_0/\mathrm{}\kappa +n\lambda )}{\kappa }})_c`$
Now suppose a number projective measurement is made on the auxiliary condensate, projecting it to $`|m`$. The normalized conditioned state of the central condensate is

$`|\mathrm{\Phi }_c={\displaystyle \frac{{}_{a}{}^{}m|\mathrm{\Phi }_{ac}^{}}{\sqrt{Tr\{|_am|\mathrm{\Phi }_{ac}|^2\}}}}`$ (38)
which for $`\lambda =(E_0/\mathrm{}\kappa )=\kappa /2`$ turns out to be

$`|\mathrm{\Phi }_c`$ $`=`$ $`{\displaystyle d^2\beta b_\beta (A|\beta e^{im\pi }+B|\beta e^{im\pi })_c},`$ (39)
and now, if $`m`$ is even the central condensate mode is left in the required state, if it is not the experiment is to be repeated. Thus the number of atoms in the auxiliary condensate mode determines the parity of the central condensate mode. In Fig. 2 we plot the efficiency of this proccess by summing over all the even number probabilities

$$P_{even}=\underset{m=even}{}P(m)=\underset{m=even}{}Tr\{|_am|\mathrm{\Phi }_{ac}|^2\},$$
(40)
for the auxiliary mode prepared in three distinct states. The dashed line is for the auxiliary mode prepared in a number state with average number of atoms $`n`$. As it is expected in this case the probability of success is constantly 0.5 indepently of the atom number. The dotted line is for a coherent state. We see that when the auxiliary mode is prepared in this state the probability of success also attains the limit of 0.5, unless the auxiliary mode has less than one atom in average as it attains the vaccum state. However the state that presents the best efficiency is the squeezed vacuum state  as depicted by the solid line in Fig. 2, by varying the squeezing parameter $`r`$ in the top axis of the plot. We see that for $`r1.425`$ the probability of success of the event is higher than 0.5. This is probably the best situation for parity control.

### B Displacement

Now, to transform the state (32) into (31), we choose a virtual displacement on the central field, defined as follows. Consider an atomic beam acting as a displacement operator over the central mode

$`D_\delta |\mathrm{\Phi }_c`$ $`=`$ $`{\displaystyle }d^2\beta b_\beta \{\mathrm{cos}[Im(\delta ^{}\beta )](A|\beta +\delta B|\beta +\delta )`$ (42)
$`+i\mathrm{sin}[Im(\delta \beta )](A|\beta +\delta +B|\beta +\delta )\}.`$
Suppose $`\delta \mathrm{}`$ and $`|\delta ||\beta |`$. If $`\delta =(l+1/2)\pi /Im(\beta )`$, for $`l=0,1,2\mathrm{}`$, the central mode state is left in the required state (but a global phase of no importance). The displacement operation can also be directly given by quadrature ($`X=b+b^{}`$) measurement through the homodyne detection described in Sec. IV. As the parameter $`|\delta ||\beta |`$ is known to be very small, the displacement operator, is given by

$$D_\delta =e^{i|\delta |X}1+i\delta X;[D_\delta ,X]=0.$$
(43)
Knowing $`\delta `$, the measurement of $`X`$ gives the required displacement. The efficiency of this proccess is dependent upon the full knowledge of the above constants, and thus upon experimental mastering.

Rotations such as that proposed above can be realised by virtual displacement. Obviously, the state (34) can be transformed into (31) by sequential applications of the displacement and parity operations. With this procedure, the teleportation protocol is complete. Notice however that none of the above operations are unitary, being dependent on selective measurements, and thus irreversible. A figure of merit of the whole proccess can be given by adding the probability of succes of each event and dividing by the number of equiprobable events,

$$P=\left(1+P_{even}+P_D+P_{even}P_D\right)/4,$$
(44)
where $`P_D`$, the probability of success of the displacement operation, is determined by the experimental control. In the right axis of Fig. 2 we compare $`P`$ for the three states considered above for the auxiliary mode, fixing $`P_D=1`$. Again we observe that for $`r1.425`$ the squeezed vacuum is the best state for the auxiliary mode to be prepared. In Fig. 3 we analyse $`P`$ for the squeezed vacuum by varying both $`r`$ and $`P_D`$. The efficiency decreases considerably (less than $`0.3`$) when both, the squeezing parameter is high and the displacement process efficiency is low. But, when the squeezing parameter is fixed to 0, or alternatively $`P_D`$ is fixed to 1, either situations have efficiency higher than $`0.5`$ attaining the ideality for $`P_D=1`$ and $`r=0`$, the vacuum state. This situation corresponds to an empty mode. Thus if the auxiliary mode can be initially prepared in a vacuum state the protocol has high probability of success.

## VI physical implementation

A specific physical implementation of the time dependent potential can be designed, similarly to what is presented in  on optical lattices. Let us consider a picture of an atom driven on a $`|J=1/2|J=3/2`$ transition by a one-dimensional optical lattice red-detuned in 1D lin-angle-lin configuration. The optical field can be written as a superposition of opposite helicity standing waves

$`𝐄_L(z)`$ $`=`$ $`\sqrt{2}E_1[e^{i\theta /2}\mathrm{cos}(k_Lz+\theta /2)𝐞_+`$ (46)
$`+e^{i\theta /2}\mathrm{cos}(k_Lz\theta /2)𝐞_{}],`$
for a convenient choice of relative phase between the beams. The potential for atoms in the ground state is

$`U(z)`$ $`=`$ $`{\displaystyle \frac{2U_1}{3}}\{2[1+\mathrm{cos}\theta \mathrm{cos}(2k_Lz)]I`$ (48)
$`+[\mathrm{sin}\theta \mathrm{sin}(2k_Lz)]\sigma _z\}{\displaystyle \frac{\mathrm{}}{2}}\gamma 𝐁\sigma ,`$
where $`U_1`$ is the light shift produced by a single beam of amplitude $`E_1`$, while $`𝐁`$ is a magnetic field and $`\{I,\sigma _i\}`$ are the identity and Pauli spin operators in the ground-state manifold. By varying $`\theta `$ the peak-peak modulation depth of the potential and the distance between the $`|m=1/2`$ and $`|m=1/2`$ potential wells are changed by

$$U_p=\frac{4}{3}U_1\sqrt{3\mathrm{cos}^2\theta +1},k_l\mathrm{\Delta }z=\mathrm{tan}^1\left(\frac{\mathrm{tan}\theta }{2}\right),$$
(49)
respectively, while changing the longitudinal component of $`𝐁`$ shifts the minima of these wells. The transverse component of $`𝐁`$ breaks the degeneracy of the bipotential at positions of linearly polarized light. By varying $`\theta `$ appropriately one can design adiabatically time dependent potentials . By loading the optical potential with the condensate at $`\theta =\pi /2`$, and changing it from $`\pi /2`$ to $`5\pi /2`$ adiabatically, we obtain the sequential controlled collision we required, as pictured in Fig. 4. The three condensate modes should be loaded in neighbouring wells in order to produce the required protocol.

An alternative approach can be designed much in the same way of the controlled entanglement generation in atomic ensembles of Jaksch et al. , with no need of a magnetic field, but where different internal spin components are selected to promote the minima of the potential to move in relation to each other. Adiabaticity here, is then related to the speed of the the minima displacement. This alternative has an advantage in relation to the first one, which is that the condensate modes do not need to be loaded in neighbouring wells, as longer as the condensate modes 1 and 3 have the same spin component, the mode 2 with an opposite component, can be moved in relation to the other two modes effecting the required operations, despite their (mode 1 and 3) relative position.

## VII Discussion

In summary, we have speculated about teleportation of a BEC state between modes in a three-well potential (over short distances). Hard-sphere cross-collision is used as a resource to generate entanglement between modes. For the final state measurement, the protocol employs a homodyne detection of the BEC quadrature, in which Josephson coupling of the condensate mode to an auxiliary mode in a neighbouring well plays the role of atomic beam-splitter.

The teleportation protocol itself was constituted of three stages, (i) entanglement formation; (ii) measurement; and (iii) receiver operations. In (i) we limited ourselves to the situation of adiabatic time varying potential ($`\frac{dV}{dt}(E^{}E_0)/\mathrm{}`$) to avoid excitations to other states. In (ii) the homodyne detection of the BEC quadrature guarantees that for ($`N\mathrm{\Omega }/\kappa `$), a well-diluted atomic gas and short time of interaction, a condensate mode quadrature can be determined by the difference of atoms between the central mode and an auxiliary mode. In (iii) parity and displacement operations were proposed to fix the teleported state at the receiver station. Since those operations are dependent upon selective measurements, the whole teleportation protocol is limited by the operations efficiencies. The efficiency of the parity operation is given by the even-parity of number of atoms present in the auxiliary state, which can be higher than 0.5 once the auxiliary condensate can be prepared approximately to a vacuum squeezed state with squeezing parameter $`r1.425`$. Together with the virtual displacement operations, those parity operations have their efficiency limited by the full knowledge of the experimental parameters. Such is the case for the squeezing parameter, $`r`$, and the virtual displacement $`\delta =(l+1/2)\pi /Im(\beta )`$, for $`l=0,1,2\mathrm{}`$. The efficiency of the whole operation vary from 0.3 to 1 for auxiliary mode prepared in the squeezed vacuum, by varying $`r`$ and the probability of success of the displacement operation.

We stress that the proposed protocol is idealised in that the above measurement processes will be very sensitive to the presence of dissipation, such as collisions with non-condensate atoms . We expect that, together with the requirement of a small number of particles at the measurement stage ($`N\mathrm{\Omega }/\kappa `$), decoherence effects impose a severe limitation to the “size” of the BEC to be teleported. The competition between Josephson coupling and cross-collision could also represent a limitation for the above protocol. We know that if the modes are prepared in coherent states this does not represent any problem, once the Josephson coupling just adds a phase to the state. However for other states this may not be so. Thus the competition of those evolutions need to be further investigated.

As a last remark, throughout the paper we have assumed a coherent basis representation for the condensate mode states, which are better adapted to our discussion. Although not specifically adressing to coherent states for the modes, those states could also be considered, if each mode were actually entangled with another condensate to form an entangled coherent state. The scheme could then proceed as described and would effectively be entanglement swapping.

Despite the idealisation, we hope that the above considerations bring some contribution to the realisation of matter field state teleportation in the near future.

###### Acknowledgements.

The author thanks G. J. Milburn for encouraging comments, and M. H. Y. Moussa and G. A. Prataviera for enlightening discussions. This work is supported by FAPESP-Brazil, under projects 01/00530-2 and 00/15084-5.