Giovanni Mazzarella

Marco Moratti

Luca Salasnich

Flavio Toigo

Nonlinear quantum model for atomic Josephson junctions with one and two bosonic species

Keywords: numbers: 03.75.Lm, 03.75.Mn, 03.75.Kk

We study atomic Josephson junctions (AJJs) with one and two bosonic species confined by a double-well potential. Proceeding from the second quantized Hamiltonian, we show that it is possible to describe the zero-temperature AJJs microscopic dynamics by means of extended Bose-Hubbard (EBH) models, which include usually-neglected nonlinear terms. Within the mean-field approximation, the Heisenberg equations derived from such two-mode models provide a description of AJJs macroscopic dynamics in terms of ordinary differential equations (ODEs). We discuss the possibility to distinguish the Rabi, Josephson, and Fock regimes, in terms of the macroscopic parameters which appear in the EBH Hamiltonians and, then, in the ODEs. We compare the predictions for the relative populations of the Bose gases atoms in the two wells obtained from the numerical solutions of the two-mode ODEs, with those deriving from the direct numerical integration of the Gross-Pitaevskii equations (GPEs). Our investigations shows that the nonlinear terms of the ODEs are crucial to achieve a good agreement between ODEs and GPEs approaches, and in particular to give quantitative predictions of the self-trapping regime.

Introduction

The prediction [1] of Bose-Einstein condensation (BEC) and the experimental achievement of BEC [2] has played a crucial role for theoretical and experimental developments in the physics of ultracold atoms. The study of the atomic counterpart [3,4,5,6,7] of the Josephson effect which occurs in superconductor-oxide-superconductor junctions [START_REF] Barone | Physics and Applications of the Josephson effect[END_REF] -which is an example of macroscopic quantum coherence -represents one of these developments. Albiez et al. [START_REF] Albiez | [END_REF] have provided the first experimental realization of the atomic Josephson junction (AJJ) previously analyzed theoretically in a certain number of papers [3,4,5,6,7]. In 2007 Gati et al. [10] reviewed the experiment by Albiez et al. [START_REF] Albiez | [END_REF] and compared the experimental data with the predictions of a manybody two-mode model [11] and a mean-field description. In the above references the analysis of AJJs physics is carried out in the presence of a single bosonic component. The possibility to tune intra-and inter-species interactions [12,13] by means of the Feshbach resonance technique makes possible to study of AJJs with two bosonic species trapped together by double-well potentials and to use BECs mixtures as powerful instruments to investigate quantum coherence and nonlinear phenomena, with particular attention to the existence of self-trapped modes and intrinsically localized states.

In the superfluid regime the dynamics of the relative populations and relative phases of the Bose condensed atoms can be described by Josepshon's two-mode equations, which are ordinary differential equations (ODEs), see for example Refs. [5,14,15,16,17]. This description is achieved in the presence of a confining double-well potential, with a single bosonic component [5] and also with bosonic mixtures [14,15,16,17,18]. One of the most interesting aspects of AJJs analysis is to compare the predictions deriving from the ODEs with the ones obtained from the Gross-Pitaevskii equations (GPEs). For single component condensates, Salasnich et al. [6] have shown that a good agreement exists between the results obtained from the GPE and those of the ODEs. Similar agreement was obtained in [7] for AJJs realized with weakly interacting solitons localized in two adjacent wells of an optical lattice. However, the situation may be quite different for multicomponent condensates, due to the interplay of intra-and inter-species interactions which enlarges the number of achievable states (for instance, mixed symmetry states can exist only in the presence of the inter-species interaction) as well as their stability, giving to the system many more dynamical possibilities. Recently it was shown that for the two components case the integration of the ODEs allows to predict the analogous of the macroscopic quantum self-trapping phenomenon observed in AJJs with one bosonic component [15,17]. This phenomenon has been discussed for a two-components nonlinear Schrödinger model with a double-well potential by Wang and co-workers [19]. More recently, a comparison between the reduced ODEs system and the full GPE dynamics was performed, showing that, for various conditions, a good agreement exists between the two kinds of predictions [17].

The aim of the present work is to analyze how the accuracy of the two-mode approximation can be improved by taking into account the usually-neglected nonlinear terms. These terms derive from the overlaps between wave functions localized in different wells. Both for single component and for two components AJJs -introduced in the second section -we proceed from a full second quantized description of the system. In Sec. III we describe the system by the extended Bose-Hubbard (EBH) Hamiltonian. In the single component case, the EBH Hamiltonian is the two-sites restriction of the Hamiltonian considered in Refs. [20,21] to analyze bosons loaded in one dimensional optical lattices. In the two species case, the EBH Hamiltonian is the extended version of the one considered by Kuklov and Svistunov in Ref. [22] to study the counterflow superfluidity of two-species ultracold atoms. We note that the study of the two components bosonic system proceeding from a pure quantum approach is a subject of wide interest. In fact, this topic is dealt with in certain regions of the phase space in Ref. [23] and in the case of hardcore bosons as discussed in Ref. [24].

The EBH Hamiltonian sustains the dynamics of the single-particle operators via the Heisenberg equations of motion [25,26]. By performing the mean-field approximation on the single-particle operators of each component, the improved ODEs are achieved. In the third section we also discuss how it is possible to distinguish the Rabi regime, the Josephson regime, and the Fock regime. This analysis is carried out in terms of the macroscopic parameters involved in the EBH Hamiltonians and, then, at the right hand sides of the improved ODEs as discussed for single AJJs in Ref. [27]. In Sec. IV we write down the GPEs for the one and the two components AJJs. Here we compare the results obtained by numerically integrating the GPEs with the predictions obtained by numerically solving the improved ODEs. Moreover, in the fourth section we plot the phase-plane portraits of the dynamical variables fractional imbalance-relative phase. Finally, in Sec. V we draw our conclusions.

The system

We consider two interacting dilute and ultracold Bose gases denoted below by 1 and 2. We suppose that the two gases are confined in a double-well trap produced, for example, by a far off-resonance laser barrier that separates each trapped condensate in two parts, L (left) and R (right). We assume, moreover, that the two condensates interact with each other and that the trapping potential V trap (r) for both components is taken to be the superposition of a strong harmonic confinement in the radial (x-y) plane and of a double-well (DW) potential in the axial (z) direction. We model the trapping potential as:

V trap (r) = m i ω 2 i 2 (x 2 + y 2 ) + V DW (z) , (1)
where m i is the mass of the ith component. For simplicity we take ω 1 = ω 2 ≡ ω. For symmetric configurations in the z direction, we take -for the ith species -the double-well in Eq. (1) as

V DW (z) = V L (z) + V R (z) V L (z) = -V 0 Sech 2 ( z + z 0 b ) V R (z) = -V 0 Sech 2 ( z -z 0 b ) V 0 = ω i [1 + Sech 2 ( 2z 0 b )] -1 , (2)
that is the combination of two Pöschl-Teller (PT) potentials, V L (z) and V R (z), centered at the points -z 0 and z 0 , and separated by a potential barrier which may be changed by varying b (see Fig. 1). We use PT potentials only for the benefit of improving accuracy in our numerical GPEs calculations (see the fourth section), taking advantage of the integrability of the underlying linear system. We remark, however, that our results apply to a generic double-well potential. Eigenvalues and eigenfuctions of the Pöschl-Teller potential for a single well are known analytically. The wave functions of the ground state of V α (z) (α = L, R), centered around -z 0 (+z 0 ) are [START_REF] Landau | Quantum Mechanics: Non-Relativistic Theory[END_REF] :

φ (α,i,P T ) (z) = A[1 -T anh 2 ( z ± z 0 b )] B i /2 B i = - 1 2 + 2m i V 0 b 2 2 + 1 4 . (3)
The constant A, in Eq. ( 3), ensures the normalization of the wave function in each well.

The second quantization Hamiltonian

To describe our system at zero-temperature, we proceed from the second quantized Hamiltonian, which reads

Ĥ = i=1,2 d 3 r Ψ † i (r) - 2 2m i ∇ 2 + V trap (r) Ψi (r) + i=1,2 g i 2 d 3 r Ψ † i (r) Ψ † i (r) Ψi (r) Ψi (r) + g 12 d 3 r Ψ † 1 (r) Ψ † 2 (r) Ψ2 (r) Ψ1 (r) , (4)
where V trap (r) is the potential (1). The coupling constants g i and g 12 are the intra-and inter-species atom-atom interaction strengths, respectively. These constants are given by

g i = 4π 2 a i m i , (5)
g 12 = 2π 2 a 12 m r , (6)
where the reduced mass m r is equal to m 1 m 2 /(m 1 + m 2 ). Eqs. ( 5) and ( 6) relate the two coupling constants to the respective s-wave scattering lengths, a i and a 12 .

In the following, we shall consider both g i and g 12 as free parameters, due to the possibility of changing the s-wave scattering lengths a i and a 12 by the technique of Feshbach resonances. In the following, we will neglect the mass difference between the two bosonic components of the mixture, as for example in Ref. [13], and assume that m 1 = m 2 ≡ m. In Eq. ( 4), the field Ψi (r) ( Ψ † i (r)) destroys (creates) a boson of the ith species at the point r, and obeys the usual bosonic commutation relations. We expand the field operator Ψi (r) in terms of operators âα,i (â † α,i ) -destroying (creating) a boson of the ith species in the well α = L, R -according to:

Ψi (r) = α=L,R Φ α,i (r)â α,i , (7)
where â's and â † 's satisfy the usual boson commutation relations and the functions Φ α,i form an orthonormal set. Due to the form (1) of the trapping potential, Φ α,i (r) can be decomposed as

Φ α,i (r) = w i (x)w i (y)φ α,i (z) , (8)
where w i (x) and w i (y) are the ground state wave functions of the harmonic oscillator potentials m i ω 2 i x 2 /2 and m i ω 2 i y 2 /2, respectively. The functions φ L,i (z) and φ R,i (z) at right hand side of Eq. ( 8) are two functions well localized in the left and right well, respectively. These functions are real and orthonormal. The functions φ L,i (z) and φ R,i (z) can be determined following the same perturbative approach as in Ref. [17].

Under the same conditions, these functions may be written in terms of the φ (L,i,P T ) (z) and φ (R,i,P T ) (z) of Eq. ( 3) as:

φ L,i (z) = 1 2 ( 1 √ 1 + s + 1 √ 1 -s )φ (L,i,P T ) (z) + ( 1 √ 1 + s - 1 √ 1 -s )φ (R,i,P T ) (z) φ R,i (z) = 1 2 ( 1 √ 1 + s - 1 √ 1 -s )φ (L,i,P T ) (z) + ( 1 √ 1 + s + 1 √ 1 -s )φ (R,i,P T ) (z) , (9)
where s = +∞ -∞ dz φ (L,i,P T ) (z)φ (R,i,P T ) (z).

AJJs with a single bosonic species

Let us start our analysis by considering the presence of a single bosonic component.

In this case, the inter-species coupling constant ( 6) is equal to zero. We use the field operator expansion [START_REF] Barone | Physics and Applications of the Josephson effect[END_REF] in the second quantized Hamiltonian (4). The AJJs microscopic dynamics is controlled by the EBH Hamiltonian [20,25,26]. The EBH model, by omitting the species index i, is described by the Hamiltonian

ĤEBH = E 0 L â † L âL + E 0 R â † R âR + U L 2 â † L â † L âL âL + U R 2 â † R â † R âR âR -K(â † L âR + â † R âL ) + K c (â † L nL âR + â † L nR âR + â † R nL âL + â † R nR âL ) + V â † L â † R âL âR + K p (â † L â † L âR âR + â † R â † R âL âL ) . (10)
Here nα = â † α âα is the number of particles in the αth well. E 0 α are the energies of the two wells, U α > 0 are the boson-boson repulsive interaction amplitudes, and K is the tunnel matrix element, which is the Rabi oscillation energy in the case of a model with U α equal to zero. The parameter K c is the induced collisionally hopping amplitude, V is the density-density bosonic interaction amplitude, and K p describes the pair bosonic hopping [20]. By using the decomposition (8) and the explicit form of w(x) and w(y), the macroscopic parameters (10) may be shown to be related to the intra-species coupling constant (5) and to the other microscopic parameters (the mass and the frequency of the harmonic trap) by the formulas

E 0 α = dz 2 2m ( dφ α dz ) 2 + (V DW + 2 2ma 2 ⊥ + mω 2 a 2 ⊥ 2 )(φ α ) 2 U α = g +∞ -∞ dz (φ α (z)) 4 K = -dz 2 2m dφ L dz dφ R dz + V DW (z)φ L φ R K c = g +∞ -∞ dz (φ α (z)) 3 φ β (z) V = 2g +∞ -∞ dz (φ α (z)) 2 (φ β (z)) 2 K p = V 4 , (11)
where g = g 2πa 2

⊥

. We observe that the first two lines of Hamiltonian (10) involve only the overlaps between φ α 's localized in the same well, see [5]. The third and fourth lines of the Hamiltonian (10) include also the overlaps between φ α 's localized in different wells, see [27]. Proceeding from the Hamiltonian (10), we write down criteria to individuate different oscillations regimes sustained by the AJJs dynamics. To this end, as discussed in Refs. [4,27], we express the Hamiltonian (10) in terms of the following operators:

Ĵx = 1 2 (â † L âL -â † R âR ) Ĵy = i 2 (â † L âR -â † R âL ) Ĵz = 1 2 (â † L âR + â † R âL ) , (12)
and the SU(2) algebra invariant Ĵ2 = ( N /2)( N/2 + 1), with N being equal to nL + nR [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF].

We assume that the two potential wells are symmetric,

E 0 L = E 0 R ≡ E and U L = U R ≡ U .
Neglecting constant terms, and using the fact that N 1, we get

Ĥ = (U -V ) Ĵ2 x -2 (K -K c N ) Ĵz + V Ĵ2 z . (13)
Here, if the condition (U -V )

V is verified, we can consider only the terms in Ĵ2

x and Ĵz . Then by defining the parameter R as

R = (U -V ) N (K -K c N ) ( 14
)
we are able (see Ref. [4]) to distinguish the three following regimes

• Rabi: R 1;

• Josephson: 1 R N 2 ;

• Fock: N 2 R.

In the Rabi regime the bosons are in a coherent state and oscillate with a frequency given simply by the energy difference between the ground state and the first excited state associated to the double-well potential. In the Josephson regime the bosons are in a coherent state and oscillate with a frequency which depends on the parameters U , K c , and V . Moreover, if the interaction strength is sufficiently large the self-trapping takes place. In the Fock regime the bosons are in a Fock state characterized by the suppression of number fluctuations. Now, we observe that the Hamiltonian (10) can be viewed as the two sites restriction of the Hamiltonian considered in Refs. [20,21] within the study of bosons loaded in one dimensional optical lattices. In particular, in Ref. [20] it is shown that within the Fock regime two regions open up. To this end we denote by ∆ the energetic gap between the Fock state with N 0 bosons per well and the Fock state with (N 0 + 1) per well [20]

∆ = 2(E + U N 0 ) . (15)
Then, when

|4∆ -(2 N 0 + 1) V | > V (16)
we have a pure Mott insulating phase (PMI), driven by the density-density on-site interaction. When

|4∆ -(2 N 0 + 1) V | < V (17)
we have a Density-Wave Mott insulating (DWMI) regime, driven by the nearestneighbors interaction [20]. Note that the DWMI phase is characterized by number fluctuations suppression as well.

At this point, we remark that we are interested in determining the fully coherent dynamical oscillations of population of the Bose condensed atoms between the left and right wells. Then, we proceed from the Heisenberg equations of motion for the model Hamiltonian (10). These equations of motion control the temporal evolution of âα . We observe that in the superfluid regime the system is in a coherent state and the following mean-field approximation [25]

âα = N α exp(iθ α ) nα = N α (18)
can be performed. The averages involved in Eq. ( 18) are evaluated with respect to the coherent state. Under the assumption of symmetric wells and by inserting the meanfield approximation (18) into the aforementioned Heisenberg equations of motion, we get

ż(t) = - 2(K -K c N ) 1 -z 2 (t) sin θ(t) + V N 2 (1 -z 2 (t)) sin 2θ(t) θ(t) = U -V Nz(t) + 2(K -K c N ) z(t) cos θ(t) 1 -z 2 (t) - V N 2 z(t) cos 2θ(t) , (19)
where N = N L +N R is the total number of bosons, and z = (N L -N R )/N and θ = θ R -θ L are, respectively, the fractional imbalance and the relative phase.

AJJs with two bosonic species

In this subsection we shall consider AJJs in the presence of two interacting bosonic components. In this case both the coupling constants ( 5) and ( 6) are finite, and the two mode EBH model is described by the Hamiltonian

Ĥ = i=1,2 Ĥ(EBH,i) + Ĥ12 . ( 20
)
The Hamiltonian Ĥ(EBH,i) is the single component Hamiltonian (10) written in terms of the operators âα,i and â † α,i . The parameters E 0 α , U α , K, K c , V , K p , and the function

φ α will read E 0 α,i , U α,i , K i , K c,i , V i , K p,i
, and φ α,i , respectively. The microscopic quantities referred to a single bosonic component will be modified according to the same prescription. Under the hypothesis of symmetric wells, the coupling Hamiltonian Ĥ12 reads:

Ĥ12 = U 12 (â † L,1 â † L,2 âL,1 âL,2 + â † R,1 â † R,2 âR,1 âR,2 ) + V 12 (â † L,1 â † R,2 âL,1 âR,2 + â † L,2 â † R,1 âL,2 âR,1 ) + K p,12 (â † L,1 â † L,2 âR,2 âR,1 + â † R,1 â † R,2 âL,2 âL,1 + â † L,1 â † R,2 âL,2 âR,1 + â † R,1 â † L,2 âR,2 âL,1 ) + K c,12 (â † L,1 nL,2 âR,1 + â † R,1 nL,2 âL,1 + â † L,2 nL,1 âR,2 + â † R,2 nL,1 âL,2 + â † L,2 nR,1 âR,2 + â † R,2 nR,1 âL,2 + â † L,1 nR,2 âR,1 + â † R,1 nR,2 âL,1 ) . (21)
In Eq. ( 21), U 12 is the inter-species interaction amplitude between bosons localized in the same well, and V 12 is the inter-species interaction amplitude between bosons localized in different wells. The quantity K p,12 is the inter-species pair hopping (hopping of particle-particle or hole-hole pair made up of bosons of different species); K c,12 is the amplitude of the inter-species collisionally induced hopping. By using the decomposition (8) and the explicit form of w i (x) and w i (y), the aforementioned parameters are shown to be related to the inter-species coupling constant (6) by:

U 12 = g12 +∞ -∞ dz (φ α,i (z)) 2 (φ α,j (z)) 2 V 12 = g12 +∞ -∞ dz (φ α,i (z)) 2 (φ β,j (z)) 2 K c,12 = g12 +∞ -∞ dz (φ α,i (z)) 3 (φ β,j (z)) K p,12 = V 12 , (22)
where g12 =

g 12 π(a 2 ⊥,1 + a 2 ⊥,2 )
. Note that we are considering both the overlaps between φ α 's localized in the same well (U α,i and U 12 ) -that are the only terms taken into account in Ref. [17] -and the overlaps between φ α 's localized in different wells (V 12 , K p,12 , K c,12 ). We observe that, in general, due to the presence of the parameters (22) the identification of different oscillation regimes proceeding from the Hamiltonian (20) is not immediate as for single component AJJs. Nevertheless, under certain conditions we are able to write down criteria to select the different regimes sustained by the two components AJJs dynamics. First, let us focus on the case in which only the overlaps between φ α 's localized in the same well are considered. If certain relations exist between the intra-and the inter-species interactions amplitudes, we can recognize the two-species corresponding of the Rabi, Josephson and Fock regimes discussed in the case of single component AJJs. For each component i, we define the quantity γ i as

γ i = U i N i K i . (23)
We recognize the following "weak-coupled" Rabi, Josephson, and Fock regimes

• Rabi: γ i 1, |U 12 | U i ; • Josephson: 1 γ i N 2 i , |U 12 | ≤ U i ; • Fock: N 2 i γ i .
In the Josephson regime, even if the intra-species interaction is not strong enough to ensure self-trapping by itself, self-trapping occurs when the inter-species interaction strength exceeds a crossover value. In the Fock regime the net number of atoms in the transport is suppressed. However, with repulsive inter-species interaction the so-called counterflow survives [22]. This means that the currents of the two species are equal in absolute values and are in opposite directions. This conductive regime is named super(counter)fluid phase (SCF). As discussed in Ref. [22], the system supports the SCF phase of the two components when

U 1 + U 2 -2 U 12 1 . (24)
When the condition

U 1 + U 2 = 2 U 12 (25)
is met, a phase separation (PS) is observed in the system and the system can be viewed as composed by two totally independent Bose gases confined in the double-well potential. On a physical level, this phase separation means that one bosonic component will occupy the left well and the other the right well. If the inter-species interaction is attractive and the hypothesis N 1 = N 2 ≡ N is verified, then, when

U 1 + U 2 -2 |U 12 | 1 ( 2 6 )
a superfluid phase, in which the superfluid consists of pairs of bosons, is supported by the system. This phase is named superfluid paired phase [START_REF] Kuklov | [END_REF]. So far we have neglected the role played by the terms deriving from the overlaps between φ α 's localized in different wells. The presence of these terms makes the scenario more complicated. However, also in this situation, under certain conditions, it is possible to achieve a classification of the oscillations regimes. To this end, as discussed for the single component case, we express the Hamiltonian (20) in terms of the operators Ĵx,i , Ĵy,i , Ĵz,i defined in Eq. ( 12) and the SU(2) algebra invariant Ĵ2 i = ( Ni /2)( Ni /2 + 1), with Ni being equal to nL,i + nR,i . Since we are assuming symmetric potential wells, we can write that E 0

L,i = E 0 R,i ≡ E i , U L,i = U R,i ≡ U i .
Neglecting constant terms, and using the fact that

N i 1, we get Ĥ = (U i -V i ) Ĵ2 x,i -2 (K -K c,i N i -K c,12 N j ) Ĵz,i + V i Ĵ2 z,i + 4((U 12 -V 12 ) Ĵx,1 Ĵx,2 + V 12 Ĵz,1 Ĵz,2 ) + U 12 (n L,1 nR,2 + nR,1 nL,2 ) + V 12 (n L,1 nL,2 + nR,1 nR,2 ) . (27) Again, if (U i -V i ) V i , V 12
, and (U 12 -V 12 ) V i , V 12 , we can consider only the terms in Ĵ2

x,i , Ĵz,i , and Ĵx,1 Ĵx,2 . We will assume also that

N 1 = N 2 ≡ N , U 1 = U 2 ≡ U , V 1 = V 2 ≡ V , and K c,1 = K c,2 ≡ K c ,
and that the initial conditions are the same for both the components. In analogy to the case of a single component AJJ, we define the parameter R as

R = ((U -V ) + 4(U 12 -V 12 )) N (K -(K c + K c,12 )N ) . ( 28
)
Again, we are able to distinguish the three regimes:

• Rabi: R 1;

• Josephson: 1 R N 2 ;

• Fock:

N 2 R.
At this point, we remark that we are interested in determining the fully coherent dynamical oscillations of population of the two bosonic components between the left and right wells. Then, we proceed from the Heinseberg equations of motion for the model Hamiltonian (20). These equations of motion control the temporal evolution of âα,i . Again, by inserting the mean-field approximation valid in the superfluid regimeâα,i = N α,i exp(iθ α,i ), nα,i = N α,i -into the aforementioned Heisenberg equations of motion, one gets the coupled differential equations for the fractional imbalance

z i = (N L,i -N R,i )/N i and relative phase θ i = θ R,i -θ L,i of the two species: żi (t) = - 2(K i -K c,i N i ) 1 -z 2 i (t) sin θ i (t) + V i N i 2 (1 -z 2 i (t)) sin 2θ i (t) + 2 (V 12 1 -z 2 j (t) cos θ j (t) + K c,12 )N j 1 -z 2 i (t) sin θ i (t) θi (t) = U i -V i N i z i (t) + 2(K i -K c,i N i ) z i (t) cos θ i (t) 1 -z 2 i (t) - V i N i 2 z i (t) cos 2θ i (t) + U 12 -V 12 N j z j (t) - 2 (V 12 1 -z 2 j (t) cos θ j (t) + K c,12 )N j z i (t) cos θ i (t) 1 -z 2 i (t) . ( 29
)

Gross-Pitaevskii equations predictions: comparison with ordinary differential equations results

So far we have discussed how AJJs dynamics can be described by means of the ODEs, i.e. Eqs. ( 19) and [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF]. We know that AJJs dynamics can be analyzed, in the mean-field approximation, in terms of partial differential equations, i.e. the GPEs. This description can be achieved proceeding from the Heisenberg motion equations for the field operators Ψi (r, t), (i = 1, 2), associated to the Hamiltonian (4), that is

i ∂ t Ψi = [ Ψi , Ĥ] . (30)
The average -denoted by ... -of both sides of Eq. ( 30) evaluated with respect to the coherent state, provides the two coupled GPEs

i ∂Ψ i ∂t = - 2 2m i ∇ 2 Ψ i + [V trap (r) + g i |Ψ i | 2 + g ij |Ψ j | 2 ]Ψ i . (31)
The macroscopic wave functions Ψ i (r, t) = Ψi (r, t) of interacting BECs in the trapping potential V trap (r) at zero-temperature satisfy Eq. ( 31). The wave function Ψ i (r, t) is subject to the normalization condition

d 3 r |Ψ i (r, t)| 2 = N i . (32)
We are interested to study the dynamical oscillations of the populations of each condensate between the left and right wells when the the barrier is large enough so that the link is weak. To exploit the strong harmonic confinement in the (x-y) plane and get the effective one-dimensional (1D) equations describing the dynamics in the z directions, we write the Lagrangian associated to the GPE equations in (31)

L = d 3 r i=1,2 Ψi (i ∂ ∂t + 2 2m i ∇ 2 )Ψ i -V trap (r)|Ψ i | 2 - g i 2 |Ψ i | 4 -g ij |Ψ i | 2 |Ψ j | 2 , (33)
where Ψi denotes the complex conjugate of Ψ i , and i = j; then, by following the decomposition (8) and the Gaussian approximation for the radial part of wave function, we adopt the ansatz

Ψ i (x, y, z, t) = 1 √ πa ⊥,i exp - x 2 + y 2 2a 2 ⊥,i f i (z, t) , (34)
where the field

f i (z, t) obey to +∞ -∞ dz|f i (z)| 2 = N i ,
so that the normalization condition given by Eq. ( 32) is satisfied. Note that the Gaussian ansatz with the transverse width simply given by a ⊥,i is reliable under very strong transverse confinements, namely when gi |f i | 2 2 ω i [START_REF] Salasnich | [END_REF]. By inserting the ansatz (34) in Eq. ( 33) and performing the integration in the radial plane, we obtain the effective 1D Lagrangian for the field f i (z, t). Such an effective 1D Lagrangian reads 2) are chosen to be b = 1 and z 0 = 3. The dashed line represents data from the integration of GPE (37), the continuous line represents data from the integration of ODEs (19), and the dot-dashed one represents data from the integration ODEs (19) with K c = V = 0. We have set N = 200 and K = 4.955 × 10 -3 . We have used the initial conditions z(0) = 0.6 and θ(0) = 0. In the top panels (from left to right): U = 0.05 K, K c = -1.842 × 10 -6 , V = 2.268 × 10 -7 ; U = 0.1 K, K c = -3.684 × 10 -6 , V = 4.535 × 10 -7 . In the bottom panels (from left to right): U = 0.2 K, K c = -7.368 × 10 -6 , V = 9.070 × 10 -7 ; U = 0.5 K, K c = -1.842 × 10 -5 , V = 2.268 × 10 -6 . Time is measured in units of ω -1 and energies are measured in units of ω.

L = dz i=1,2 fi (i ∂ ∂t + 2 2m i ∂ 2 ∂z 2 )f i -( i + V DW (z))|f i | 2 - gi 2 |f i | 4 -gij |f i | 2 |f j | 2 , (35
where i is given by i =

2 2m i a 2 ⊥,i + m i ω 2 i a 2 ⊥,i 2
. By varying L with respect to fi , we obtain the 1D GPE for the field 19) with U = 0.05 K and K c = V = 0; the continuous line represents data from the ODE (19) with U = 0.05 K, K c = -1.842 × 10 -6 and V = 2.268 × 10 -7 . The right panel shows the phase diagram for the self-trapping. In this panel the dashed line represents data from the ODE (19) with U = 0.2 K and K c = V = 0; the continuous line represents data from the ODE (19) with U = 0.2 K, K c = -7.368×10 -6 and V = 9.070 × 10 -7 . Initial conditions are the same as in Fig. 2. Time is measured in units of ω -1 and energies are measured in units of ω.

f i i ∂f i ∂t = - 2 2m i ∂ 2 f i ∂z 2 + [ i + V DW (z) + gi |f i | 2 + gij |f j | 2 ]f i . (36)
In the presence of a single bosonic component, g 12 = 0; then, the two coupled 1D GPEs Eq. (36), omitting the species index i, reduce to

i ∂f ∂t = - 2 2m ∂ 2 f ∂z 2 + [ + V DW (z) + g|f | 2 ]f . (37)
Now, we observe that it is possible to write the fields f i , (i = 1, 2), by using the two-mode approximation as done, for example, in Ref. [17]

f i (z, t) = ψ L,i (t)φ L,i (z) + ψ R,i (t)φ R,i (z) ψ α,i (t) = N α,i (t) exp(iθ α,i (t)) , (38)
with φ α,i (z) constructed as discussed in Sec. III, see Eq. ( 9). Then, one takes into account the overlaps both between φ α 's localized in the same well and between φ α 's localized in different wells. By following the same path as in Ref. [17], when the interspecies coupling constant g 12 is finite, it is possible to recover Eqs. [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF] for binary AJJs, while for g 12 equal to zero one gets back the Eqs. (19) for single component AJJs . At this point -both for single component and for two components AJJs -we may compare the predictions of the ODEs, Eqs. ( 19) and [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF], and those of the GPEs, Eqs. (37) and (36). The results of this analysis are reported in Fig. 2 for the single component case, and in Figs. 4,6 for the two components case. In obtaining Fig. 2 we have fixed the parameters b and z 0 of the double-well potential (2). Then, by using the functions (9) into the third of Eqs. (11), we have obtained the tunneling amplitude 1 0 1 2 1.0 0.5 0.0 0.5 1.0

Θ 1 t Π z 1 t 2 1 0 1 2 1.0 0.5 0.0 0.5 1.0 Θ 2 t Π z 2 t Figure 5.
Phase diagrams of the fractional imbalance z i (t) vs. macroscopic phase θ i (t) of the two bosonic species. The parameters of the double-well potential (2) are the same as in Fig. 4. In both the panels we have set

N 1 = 200, N 2 = 100, K 1 = K 2 ≡ K = 4.955 × 10 -3 , U 1 = U 2 ≡ U = 0.1K.
In both the panels, the dashed line represents data from the ODEs (29) with U 12 = -U/2 and K c,i = V i = K c,12 = V 12 = 0, the continuous line represents data from the ODEs with

U 12 = -U/2, K c,1 = K c,2 ≡ K c = -3.684 × 10 -6 , V 1 = V 2 ≡ V = 2.268 × 10 -7 , K c,12 = -K c /2, V 12
= -V/4. Initial conditions are the same as in Fig. 4. Time is measured in units of (ω 1 ) -1 = (ω 2 ) -1 ≡ ω -1 and energies are measured in units of ω. Fractional imbalance z i (t) of the two bosonic species vs. time. In the double-well potential (2) we set b = 1 and z 0 = 3. In this figure the dashed line represents data from the integration of GPEs (36), the continuous line represents data from the integration of ODEs [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF], and the dot-dashed line represents data from the integration ODEs (29) with K c,i = V i = K c,12 = V 12 = 0 (i.e. the ODEs analyzed in Ref. [17]). We have fixed N 1 = 200 and N 2 = 100. Moreover,

K 1 = K 2 ≡ K = 4.955 × 10 -3 , U 1 = U 2 ≡ U = 0.1 K, K c,1 = K c,2 ≡ K c = -3.684 × 10 -6 , V 1 = V 2 ≡ V = 2.268×10 -7 , U 12 = -2 U , K c,12 = -2 K c V 12 = -V .
Initial conditions are the same as in Fig. 4. Time is measured in units of (ω 1 ) -1 = (ω 2 ) -1 ≡ ω -1 and energies are measured in units of ω.

intra-species interaction is not strong enough to induce oscillations of z i around a non zero time-averaged value. Nevertheless, if the inter-species interaction is sufficiently large, oscillations of z i around z i (t) = 0 are observed. We have reported this kind of behavior for both the components in Fig. 6. From this figure we can see that, especially in the case of large inter-species interaction, the role played by the parameters describing the overlaps between φ α 's localized in different wells becomes essential to improve the agreement between the ODEs and the GPEs predictions. Moreover, in Fig. 3 -single component AJJs -and in Fig. 5 and Fig. 7 -two components AJJs -we show the phase-plain portraits of the dynamical variables z i and θ i for different values of the macroscopic parameters ( 11) and ( 22) (see Figs. 3,5, 7 for the details). These figures show the comparison between the trajectories in the phases space obtained by integrating the ODEs in the absence of the usually-neglected nonlinear terms (dashed lines), and the trajectories obtained from the improved version of ODEs (continuous lines). In particular, the left and the right panels of Fig. 3 show the phases space trajectories for the Josephson and self-trapping regimes, respectively, for single component AJJs. For two components AJJs, in Fig. 5 we have plotted the phases space trajectories for the Josephson regime, and in Fig. 7 we have plotted the phases space trajectories when the system is self-trapped. From Figs. 3, 5, 7 we can see that the trajectories predicted when the ODEs Eqs. ( 19) and ( 29) are solved in the absence of the usually-neglected nonlinear terms are sufficiently close to those predicted when these ODEs are solved in the presence of the aforementioned terms. Then, the dynamical evolution predicted by the standard ODEs reveals to have a good degree of reliability.

Conclusions

We have analyzed atomic Josephson junctions for a single Bose gas and for binary mixtures of bosons in a double-well potential along the axial direction and a strong harmonic confinement in the transverse directions. We have shown that for both the cases the Hamiltonian belongs to the extended Bose-Hubbard model and besides the density-density interaction it contains the pair hopping and collisionally induced hopping terms. These terms derive from the overlaps between wave functions localized in different potential wells. We started from these Hamiltonian models and established connections with spin Hamiltonians. Proceeding from these, we have discussed the possibility to discriminate, under certain conditions, different dynamical regimes sustained by the bosonic junctions. From the mean field analysis of the equations of motion for the single-particle operators involved in the extended Bose-Hubbard Hamiltonians, we have obtained the ordinary differential equations that control the macroscopic dynamics of the atomic Josephson junctions. Within the analysis of the atomic Josephson junctions macroscopic dynamics we have plotted the phase-plane portraits of the dynamical variables (fractional imbalance-relative phase) showing that the inclusion of the aforementioned collisionally induced hopping and pair hopping terms are crucial to get good agreement between the dynamics of the Josephson model described by ordinary differential equations and the one of the time dependent Gross-Pitaevskii equations, especially when the atom-atom interaction is strong.

Finally, it is important to remark that the obtained results are of general validity also for more confining (e.g. not saturating to zero at large distances) double-well potentials. Nevertheless, it is possible to design a model of pair hopping and collisionally induced hopping for bosonic atoms that is physically meaningful when optical lattices play the role of confining potentials. Physical effects related to pair hopping and collisionally induced hopping should be observable in generalizations of current experiments to detect the superfluid and insulating phases [33].

Figure 1 .

 1 Figure 1. The double-well potential (2) as a function of z for z 0 = 3 and different values of b. The dot-dashed line corresponds to b = 0.7, the continuous line corresponds to b = 1, and the dashed line corresponds to b = 1.3. Lengths are measured in units of a ⊥,i = m i ω and energies in units of ω .

Figure 2 .

 2 Figure 2. Fractional imbalance z(t) vs. time for single component atomic Josephson junctions. The parameters of the double-well potential (2) are chosen to be b = 1 and z 0 = 3. The dashed line represents data from the integration of GPE (37), the continuous line represents data from the integration of ODEs(19), and the dot-dashed one represents data from the integration ODEs(19) with K c = V = 0. We have set N = 200 and K = 4.955 × 10 -3 . We have used the initial conditions z(0) = 0.6 and θ(0) = 0. In the top panels (from left to right): U = 0.05 K, K c = -1.842 × 10 -6 , V = 2.268 × 10 -7 ; U = 0.1 K, K c = -3.684 × 10 -6 , V = 4.535 × 10 -7 . In the bottom panels (from left to right): U = 0.2 K, K c = -7.368 × 10 -6 , V = 9.070 × 10 -7 ; U = 0.5 K, K c = -1.842 × 10 -5 , V = 2.268 × 10 -6 . Time is measured in units of ω -1 and energies are measured in units of ω.

Figure 3 .

 3 Figure 3. Phase diagrams of the fractional imbalance z(t) vs. macroscopic phase θ(t) for single component atomic Josephson junctions. The parameters of the double-well potential (2) are the same as in Fig. 2. In both the panels we have set N = 200 and K = 4.955 × 10 -3 . Left panel: the dashed line represents data from the ODE (19) with U = 0.05 K and K c = V = 0; the continuous line represents data from the ODE (19) with U = 0.05 K, K c = -1.842 × 10 -6 and V = 2.268 × 10 -7 . The right panel shows the phase diagram for the self-trapping. In this panel the dashed line represents data from the ODE (19) with U = 0.2 K and K c = V = 0; the continuous line represents data from the ODE(19) with U = 0.2 K, K c = -7.368×10 -6 and V = 9.070 × 10 -7 . Initial conditions are the same as in Fig.2. Time is measured in units of ω -1 and energies are measured in units of ω.

Figure 6 .

 6 Figure 6.Fractional imbalance z i (t) of the two bosonic species vs. time. In the double-well potential (2) we set b = 1 and z 0 = 3. In this figure the dashed line represents data from the integration of GPEs (36), the continuous line represents data from the integration of ODEs[START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF], and the dot-dashed line represents data from the integration ODEs (29) with K c,i = V i = K c,12 = V 12 = 0 (i.e. the ODEs analyzed in Ref.[17]). We have fixed N 1 = 200 and N 2 = 100. Moreover,K 1 = K 2 ≡ K = 4.955 × 10 -3 , U 1 = U 2 ≡ U = 0.1 K, K c,1 = K c,2 ≡ K c = -3.684 × 10 -6 , V 1 = V 2 ≡ V = 2.268×10 -7 , U 12 = -2 U , K c,12 = -2 K c V 12 = -V .Initial conditions are the same as in Fig.4. Time is measured in units of (ω 1 ) -1 = (ω 2 ) -1 ≡ ω -1 and energies are measured in units of ω.

This work has been partially supported by Fondazione CARIPARO through the Project 2006: "Guided solitons in matter waves and optical waves with normal and anomalous dispersion". G. M. thanks A. B. Kuklov and B. V. Svistunov for useful comments.

Fractional imbalance z i (t) of the two bosonic species vs. time. The parameters of the double-well potential (2) are chosen to be b = 1 and z 0 = 3. Here, the dashed line represents data from the integration of GPEs (36), the continuous line represents data from the integration of ODEs [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF], and the dot-dashed line represents data from the integration ODEs (29) with K c,i = V i = K c,12 = V 12 = 0 (i.e. the ODEs analyzed in Ref. [17]). We have fixed N 1 = 200 and N 2 = 100. Moreover,

. We used the initial conditions z 1 (0) = 0.5 = -z 2 (0) and θ 1 (0) = θ 2 (0) = 0. In the top panels we set

and energies are measured in units of ω.

Phase diagrams of the fractional imbalance z i (t) vs. macroscopic phase θ i (t) of the two bosonic species for the self-trapping. The parameters of the double-well potential (2) are the same as in Fig. 4. In both the panels we have set

In both the panels, the dashed line represents data from the ODEs [START_REF]This definition of Ĵx and Ĵz follows that of Refs[END_REF] with 12 = V 12 = 0, the continuous line represents data from the ODEs with

Initial conditions are the same as in Fig. 4. Time is measured in units of (ω 1 ) -1 = (ω 2 ) -1 ≡ ω -1 and energies are measured in units of ω.

K. We have keeped fixed K and we have plotted the predictions of the ODEs (19) for z(t) in correspondence to different intra-species interactions both when K c and V are zero -dot-dashed lines -and in the presence of K c and V -continuous lines; the dashed lines represent z(t) obtained by numerically integrating the GPE (37). In Figs. 4,6 we have fixed the tunneling amplitude K i -as done previously in the single component case -and the intra-species interaction U i , and we have plotted the predictions of the ODEs (29) for z i (t) in correspondence to different inter-species interactions both when K c,i , V i , K c,12 , V 12 are all equal to zero -dot-dashed lines -and in the presence of K c,i , V i , K c,12 , V 12 -continuous lines; again, the dashed lines represent z i (t) obtained by numerically integrating the GPEs (36). In the two top panels of Fig. 2 and in all the panels of Fig. 4, we have plotted the temporal evolution of the bosonic fractional imbalances z when they oscillate around a zero time-averaged value, i.e. z(t) = 0. We see that the usually-neglected nonlinear terms play a crucial role in order to improve the agreement between the GPEs the ODEs predictions. In fact, neglecting these terms, the solutions of ODEs and GPEs diverge rather rapidly, as shown by dot-dashed lines in Fig. 2 -single component AJJs -and by dot-dashed lines in Fig. 4, for two components AJJs. The two bottom panels of Fig. 2 show the results of our analysis when the intraspecies interaction amplitude U is sufficiently large to induce oscillations of z(t) around a non zero time-averaged value, that is the self-trapping. We see that the inclusion within the description of the system of the usually-neglected nonlinear terms produces an improvement in the agreement between the ODEs and GPE predictions. In the two components case, from Fig. 4 we can see that the nonlinearity associated to the


O Fialko

K Ziegler
email: klaus.ziegler@physik.uni-augsburg.de

Noise correlations of a strongly attractive spin-1/2 Fermi gas in an optical lattice

Keywords: numbers: 68.85.-d, 05.30.Jp, 05.30.Fk

We calculate density-density correlations of an expanding gas of strongly attractive ultra-cold spin-1/2 fermions in an optical lattice. The phase diagram of the tightly bound fermion pairs exhibits a Bose-Einstein condensed state and a Mott insulating state with a single molecule per lattice site. We study the effects of quantum fluctuations on the correlations in both phases and show that they are especially important in the Bose-Einstein condensate state, leading to the appearance of singular peaks. In the Mott insulating state the correlations are characterized by sharp dips. This can be utilized in experiments to distinguish between these two phases.

Introduction

Experiments with ultra-cold atoms have opened many directions to study ground state properties of complex many-body systems, such as Bose-Einstein condensation and fermion pairing. The introduction of optical lattices has added the possibility of creating new ground states [1]. A prominent example is the formation of a Mott insulating (MI) phase for strongly interacting bosons [2] or fermions [3].

Another interesting direction was brought to the field by the unprecedented control over the interaction between atoms via a magnetic Feshbach resonance [4]. By means of the latter it became possible to turn fermions with two hyperfine internal states into bosons [5] by pairing. For weak attraction the Bardeen-Cooper-Schrieffer (BCS) theory predicts an instability towards formation of large Cooper pairs, whose size can be much larger than the separation between atoms in the Fermi gas. If the attraction is turned on stronger, then these Cooper pairs become smaller and eventually, for very strong interaction, they behave as point-like (hard-core) bosonic molecules. Although there is no sharp distinction between the BCS state and a condensate of local pairs of fermions, the cross-over from the weakly interacting to the strongly interacting Fermi gas is of substantial interest [6]. In the presence of an optical lattice the tightly bound molecules in the BEC regime may undergo a quantum phase transition from a BEC state to a MI state with one particle per lattice site. This leaves us with question how to identify the two phases by simple measurements.

A clear signature of superfluidity in a system of bosons is the presence of sharp density peaks due to phase coherence in time-of-flight experiments [2]. In the MI phase, on the other hand, the density peaks are washed out due to the absence of phase coherence (see [7]). To characterize this phase it is useful to study the noise correlations in the Fermi gas [8]. Noise correlations have been investigated for weakly interacting Fermi gases in an optical lattice. As a characteristic feature a correlation hole was observed experimentally [9].

In this paper, we study density-density correlations of an expanding gas of ultracold spin-1/2 fermions realized in an optical lattice. We discuss that the density-density correlations reveal the fermionic nature of the gas in the MI state and its bosonic nature in the Bose-Einstein condensate (BEC). We show that the effects of quantum fluctuations on the correlations are important.

Model

An effective Hamiltonian for fermions in an optical lattice near a Feshbach resonance was derived in Ref. [10]. Under typical experimental conditions, the resulting Hamiltonian can be mapped to an effective single-band model, which comprises the local dressed bosonic molecules and individual fermionic atoms separately. In this paper we consider a purely fermion model in an optical lattice with attractive interaction and with two internal states, denoted by the pseudospin projections ↑ and ↓, respectively. The number of fermions in the two components is equal. The fermionic field operator ψσ (r) is expanded in a basis of Wannier functions [11] as

ψσ (r) = i,n w n (r -R i )ĉ i,n,σ . (1)
which allows us to use the tight-binding representation. Here, ĉi,n,σ is the annihilation operator for particles in the Wannier state n at lattice site R i . For a deep optical lattice only the lowest Wannier state is taken into account with ĉi,σ ≡ ĉi,0,σ . In this case the lattice Hamiltonian in d dimensions reads

Ĥ = - t 2d σ=↑,↓ i,j ĉ † iσ ĉjσ - J 2d i,j ĉ † i↑ ĉj↑ ĉ † i↓ ĉj↓ - σ=↑,↓ i µ σ ĉ † iσ ĉiσ . (2)
Nearest-neighbor tunneling of the individual fermions is described by the parameter t. Moreover, there is also a nearest-neighbor attraction with parameter J which can be understood as a tunneling term of fermionic pairs (molecular tunneling). A similar model of a mixture of hard-core bosons and fermions in a lattice was discussed in Ref. [12]. However, in our model the hard-core bosons bi are given by pairs of fermion operators ĉi↑ ĉi↓ . This means that the fermions can transmute dynamically into hard-core bosons by local pairing.

Phase diagram

We apply the functional integral formalism developed for the Hamiltonian in Eq. ( 2) in [13,14]. The partition function of a grand-canonical ensemble of fermions with chemical potential µ at inverse temperature β can be expressed in terms of an integral over a fermionic (Grassmann) field ψ and complex molecular fields φ, χ as

Z = e -β 0 dτ L D[ψ, φ, χ] ( 3 )
with the Lagrangian

L = i,j φi v-1 i,j φ j + 1 2J i χi χ i - i,j Ĉ † i Ĝ-1 ij Ĉj , (4)
where we have used the notation Ĉ † = (ψ ↑ , ψ↓ ) and Ĉ = (ψ ↓ , ψ↑ ) T , and with the inverse fermionic Green's matrix

Ĝ-1 = -iφ -χ ∂ τ + µ + t ∂ τ -µ -t i φ + χ . ( 5
)
Here vi,j = J 2d δ |i-j|,1 . t is the nearest-neighbor tunneling matrix with elements t/2d. The integration over the complex molecular fields φ, χ can be performed in saddlepoint approximation [13,14]. For a vanishing fermionic tunneling, i.e. by assuming a situation where all fermions are paired up, the phase diagram is depicted in Fig. 1. There are three phases: the BEC of molecules with the condensed fraction n 0 = (J 2 -4µ 2 )/(4J 2 ), the MI state with one particle per site and n 0 = 0 and the empty phase. Fluctuations around the saddle point provide the low-energy excitations of the bosonic molecules. They are gapless in the BEC phase:

k = 4J 2 g k n 0 + 4g 2 k µ 2 , (6)
where g k = 1 -1/d d i=1 cos k i is the dispersion of the free Bose gas. On the other hand, the excitations of the MI state have a gap ∆ = 2µ -J > 0 :

k = ∆ + Jg k . (7)

Density-density correlations

In a typical experiment with ultra-cold atoms it is not easy to identify the spin projection. Therefore, the detected fermion density should be described by a superposition of both spin projections. For this purpose we introduce the spinindependent fermionic density n = n↑ + n↓ . Then the truncated density-density correlation function in a time-of-flight experiment is given under the assumption that the atoms expand freely after they have been released from the optical lattice as [15] C r,r ≡ n(r)n(r

) -n(r) n(r ) = M ht 2d n(k)n(k ) -n(r) n(r ) .( 8
)
Here k is related to r by k = Mr/ht with the mass of atoms M ∝ h/J. Free expansion assumes that atoms evolve independently. This can be achieved by switching off an optical lattice as well as by switching the magnetic field to values far from the strongly attractive regime suddenly before the expansion [8,16]. The above formula reveals then that the density distribution, observed after a fixed time-of-flight, is related to the momentum distribution of the particles trapped in the lattice before the expansion. It is widely used to analyze ground state properties of ultra-cold atoms trapped by an optical lattice [17]. Now we use the expansion of Eq. ( 1) to express the density-density correlation of Eq. ( 8) in terms of the coordinates of the underlying optical lattice as

C r,r = M ht 2d | w(k)| 2 | w(k )| 2 i,i ,j,j α,β=↑,↓ e ik(R i -R i )+ik (R j -R j ) ĉ † iα ĉ † jβ ĉj β ĉi α , + n(r) δ(r -r ) -n(r) n(r ) , (9)
where w is a Fourier transform of the Wannier function. Using the saddle-point approximation of the complex fields (φ, χ) we can evaluate ĉ † iα ĉ † jβ ĉj β ĉi α . Since the fermion field appears only in a quadratic form in Eq. ( 4), the integration are given by a Wick contraction of pairs of fermions:

ĉ † iα ĉ † jβ ĉj β ĉi α = ĉ † iα ĉi α ĉ † jβ ĉj β -ĉ † iα ĉj β ĉ † jβ ĉi α + ĉ † iα ĉ † jβ ĉj β ĉi α .( 10
)
The fermionic expectation values on the right-hand side are matrix elements of the fermionic Green's function Ĝ of Eq. ( 5):

ĉ † i↑ ĉj↑ = ĉ † i↓ ĉj↓ = 1 β n G 12 ij (ω n ) ( 1 1 ) and ĉ
† i↓ ĉ † j↑ = 1 β n G 11 ij (ω n ), ĉi↑ ĉj↓ = 1 β n G 22 ij (ω n ), (12)
where we have summed over Matsubara frequencies ω n . These expressions include also quantum fluctuations of the bosonic molecules when we take into account Gaussian fluctuations around the saddle-point solutions. For instance, we can integrate over these fluctuations to evaluate the correlation function

ĉ † iα ĉ † jβ ĉj β ĉi α -ĉ † iα ĉ † jβ ĉj β ĉi α . (13)
In general, the Green's function Ĝij is not diagonal in space due to the presence of the fermionic tunneling t but the decay is exponential with distance |R i -R j |. The exponential behavior can be approximated by taking only the diagonal part into account. In other words, in a strongly interacting Fermi gas the formation of bosonic molecules and their dynamics is the dominant feature, and the individual tunneling of fermions is negligible. Moreover, the Green's function depends on the fields φ and χ. In that case the density-density correlations reduces to

C r,0 = n(r) δ(r) + 2 M ht 2d | w(k)| 2 | w(0)| 2 ×   - i e ikR i ĉ † i↑ ĉi↑ 2 + i e ikR i ĉi↑ ĉi↓ 2 + S(k)   . ( 14
)
The first two terms constitute the density-density correlation function for fermions [15], whereas the third and the fourth terms account for the presence of the condensed molecules. In particular, the fourth term describes the effect of quantum fluctuations and can be expressed as a momentum distribution of the molecules [14]

S(k) = i,j e k(R i -R j ) ( b † i bj -b † i bj ) ( 1 5 ) with b † i = ĉ † i↓ ĉ † i↑ .
Phase fluctuations can destroy these terms, for instance, in the case of a MI state of the molecules. It is important to notice that fermionic terms contribute with a negative sign in Eq. ( 14), in contrast to the phase-coherent molecules, which contribute with a positive sign. This indicates a competition of the fermionic and the molecular contribution to the density-density correlation function. This provides a concept for measuring the properties of a strongly interacting Fermi gas.

Results

The expressions ĉ † i↑ ĉi↑ and ĉi↑ ĉi↓ are constant due to translational invariance:

C r,0 ∝ N 2 δ k,G -1 + N 0 N N 2 + S(k), (16)
where we have denoted 2N = 2 i ĉ † i↑ ĉi↑ as the total number of fermions and

N 0 N = | i ĉi↑ ĉi↓ | 2 ∝ N 2 |iφ + χ| 2 .
N is the number of lattice sites. The ratio n 0 = N 0 /N is the condensate fraction, i.e. the relative contribution of condensed molecules.

BEC

For momenta close to the reciprocal lattice vectors G the main contribution in the BEC comes from the term S(k), since it is singular for k ∼ 0:

C r,0 ∝ S(k) ≈ 4Jn 0 + Jg k k (1 -2g k ) . (17)
This results connects the density-density correlation function with the spectrum of the molecular condensate given in Eq. ( 6). The infrared divergence of the momentum distribution is a general property of a BEC at T = 0 [START_REF] Pitaevskii | Bose-Eisntein Condensation[END_REF]. The related 2D plot of the correlation function is shown in Fig. 2.

MI

For larger densities the phase coherence in the molecular state is destroyed and a MI state with one bosonic molecule (i.e. a pair of fermions) appears. Due to strong phase fluctuations the second term in Eq. ( 16) vanishes, and the correlation function becomes

C r,0 ∝ -N 2 δ k,G + S(k) ( 1 8 )
with the non-singular term Since the number of fermions N is large in Eq. ( 18), the contribution of S(k) is negligible in time-of-flight experiments of the MI. Thus the MI state is characterized by sharp dips in contrast to the singular peaks in the BEC, which appear on positions corresponding to the reciprocal lattice vectors G, similar to a Bragg diffraction pattern. However, in contrast to the latter, where the diffraction pattern is created by light scattering on atoms of a crystal, here the pattern is created by the atomic state itself in the expansion process of the cloud. Eq. ( 19) describes the bosonic nature of the molecules in the MI state: There are hole excitations in the MI phase (particle excitations are suppressed by the hard-core nature of the bosons) and their dynamics contributes to the interference pattern via the factor ∝ g k . A similar situation was observed experimentally in the case of bosonic Mott insulator in Ref. [START_REF] Gerbier | [END_REF]. The 2D plot of the correlation function is shown in Fig. 3.

S(k) ≈ J 2 4µ 2 (1 -2g k ). ( 19
)

Discussion

The density-density correlation function of an expanding cloud provides a clear picture of the state when it was still trapped by an optical lattice. It consists of four different contributions in Eq. ( 14), two of them are related to the fermionic nature of the atoms. The second term leads to a fermionic dip in the density-density correlation function in the MI state. A third term measures the condensate fraction in the case of condensed bosonic molecules and together with the singular fourth term leads to the sharp peaks. Even though the third term can compete with the second term, the singular term is the most relevant in the BEC phase.

The competing behavior of the fermionic dips and the condensate peaks in the density-density correlation function are a result of the anticommuting properties of the fermionic operators. This provides a simple concept to distinguish different states in a cloud of attractively interacting fermions. This behavior is rather different in a bosonic cloud, where all atoms contribute with the same sign to density-density correlation function because the bosonic operators commute [20].

In conclusion, we have studied an expanding cloud of strongly interacting spin-1/2 fermions after its release from an optical lattice. The properties are described in terms of the density-density correlation function. At lower densities a BEC is formed by the paired fermions, visible in the density-density correlation function as sharp peaks. At higher densities a Mott-insulating phase appears, characterized by dips in the densitydensity correlation function. This distinct behavior can be used in future experiments to distinguish between these two phases in a gas of strongly attractive spin-1/2 fermions in an optical lattice.

Figure 1 .

 1 Figure 1. Phase diagram for t = 0 and k B T = 0. The three phases are the BEC with a non-vanishing condensed density, the MI states with vanishing condensed density and with one molecules per lattice site and the empty phase. µ and J are in arbitrary energy units.

Figure 2 .

 2 Figure2. Density-density correlation function in units of (ht/M a) -2d of a strongly attractive spin-1/2 Fermi gas in an optical lattice with J = 1, µ = 0.4, k B T = 0, where the gas forms a BEC. The number of molecules per lattice site is 0.9. The singular behavior at the positions corresponding to the reciprocal lattice vectors is due to long range phase coherence of the condensed molecules. There is also small modulation due to the tunneling of molecules. The axes are given in units of ht/M a, where a is a lattice spacing.

Figure 3 .

 3 Figure 3. Density-density correlation function as in Fig.2but with J = 1, µ = 0.6, k B T = 0. The number of molecules per lattice site is 1. The singular behavior in Fig.2is replaced by dips here. This is due to the anticommuting nature of fermionic operators and thus reflects the fermionic nature of the MI state.


B K Sahoo
email: b.k.sahoo@rug.nl

Ab initio studies of electron correlation effects in the atomic parity violating amplitudes in Cs and Fr

Keywords: numbers: 31.15.Ar, 31.15.Dv, 31.25.Jf, 32.10.Dk atomic parity violation, ab initio method, coupled-cluster theory

HAL is

I. INTRODUCTION

One of the most challenging contemporary problems in physics is the search for possible new physics beyond the well known standard model (SM) of elementary particles [1][2][3][4]. Apart from using gigantic accelerators at high-energy scales, it is also possible to use high precision, albeit, low-energy table top atomic experiments, such as, the measurement of atomic parity violation, in combination with accurate relativistic many-body calculations of the atomic parity nonconserving transition amplitude, E 1P NC to achieve this goal [4,5]. Some of the prominent signatures of physics beyond the SM which can be inferred from these atomic experiments are: a tight limit on the mass of extra Z-bosons, precise value of the Weinberg angle, limit to the radiative corrections for the electron-nuclear weak interactions etc [4]. As some of these SM results are known to high precision, they demand similar sub-one percent accuracies in both the measurements and the atomic calculations. As the interaction Hamiltonian for the atomic parity violation (APV) due to the nuclear spin-independent (NSI) electron-nucleus interactions is proportional to Z 3 where Z being the atomic number [5], heavy atomic systems are chosen for the study of APV effects. A series of APV experiments on a number of atomic systems, including those of Cs [6] and Tl [7] have been carried out. However, the high accuracy of ∼ 0.35% has been achieved only for Cs. Furthermore, a number of ab initio calculations of APV amplitudes in Cs have been carried out using a variety of many-body approaches. Some results based on the relativistic coupled-cluster (RCC) theory are also available for Cs, however, their accuracies are somewhat uncertain, since most of them have used the sum-over-states approach which considers contributions from the core orbitals approximately and accounts for only a selected number of excited states whose contributions are dominant. In addition, the doubly excited intermediate atomic states and the normalization of the RCC wave functions are treated only approximately. We have developed a technique in the frame work of RCC theory that circumvents these drawbacks and it has been employed earlier in the calculation of the E1 P NC amplitudes in Ba + [8] and Ra + [9] in which we have demonstrated that the accuracies of < 1% and < 3%, respectively, were possible.

In this work, we employ the new RCC approach, mentioned above, to study various correlation effects in the parity violating amplitudes in Cs and Fr. We report a comparative study of their results along with those reported previously.

II. THEORY OF APV

The dominant interaction in an atom is the electromagnetic interaction which, as well known, conserves parity. However, there is a non-zero probability of the interaction between the electrons and the nucleus of an atom due to the weak force with the exchange of a Z 0 boson, as shown in Fig. 1, which violates parity. The interaction Hamiltonian between the electrons and the nucleus due to the weak interaction can have two components: one, vector-axial-vector and the other, axial-vector-vector currents. The latter depends on the nuclear spin and most of its contribution cancels out except from the odd nucleon and hence, it is relatively smaller in magnitude [5] than the former; the NSI component. In this work, we shall consider the APV effect due to the NSI component alone. The APV interaction Hamiltonian due to the NSI component is given by,

H NSI AP V = G F 2 √ 2 Q w (N )γ e 5 ρ e N (r), (2.1)
where G F is the Fermi constant, ρ e N (r) is the electron density over the nucleus and γ e 5 (= iγ e 0 γ e 1 γ e 2 γ e 3 ), is the product of the four Dirac matrices that involve electron spin, and Q w (N ) is the nuclear weak charge, which is equal to 2(Zc p + Nc d ) where c p and c n denote the electron-proton and the electron-neutron coupling constants for the atomic (Z) and neutron (N ) numbers, respectively. The values of these coupling constants predicted by the SM, in the lowest order of electroweak interaction (at the tree level), are given by

c p = 1 2 (1 -4sin 2 θ W ) ≈ 0.04, c n = - 1 2 , ( 2.2)
where θ W is the Weinberg angle and its measured value is ∼ sin 2 θ W ≈ 0.23 [10]. Substituting these values in Q w (N ), we get Q w (N ) = -N + Z(1 -sin 2 θ W ) which is proportional to N . Hence, the perturbation due to H NSI AP V is generally expressed in the scale of Q w (N )/N . Since H NSI AP V does not commute with the parity operator, its inclusion with the atomic Hamiltonian of the electromagnetic interaction, which commutes with the parity operator, mixes the opposite parity states of same angular momentum. The strength of this interaction is sufficiently weak, which justifies its consideration as a first-order perturbation.

The Dirac-Coulomb (DC) Hamiltonian used, here, in the calculation of the atomic wave functions of definite parity is given by

H DC = i cα • p i + (β -1)c 2 + V nuc (r i ) + i>j 1 r ij , (2.3)
where c is the velocity of light, α and β are the Dirac matrices (note that γ e = βα) and V nuc (r) is the nuclear potential.

The atomic wave functions (|Ψ

(0)
n ) corresponding to H DC can be considered as the unperturbed wave functions. The total wave function of a system including the first order correction due to the interaction Hamiltonian H NSI AP V is given by (1) n , (2.4) where |Ψ

|Ψ n = |Ψ (0) n + G F |Ψ
(1) n is the first order perturbed wave function of its unperturbed valence state |Ψ (0) n and G F is used as a coupling constant.

Electric dipole (E1) transitions between the states of same parity are forbidden due to the electromagnetic selection rules. However, an E1 transition between the states of mixed parity, mixed due to H NSI AP V interaction, is possible and the corresponding transition amplitude can be expressed as

E1 P NC = Ψ f |D|Ψ i Ψ f |Ψ f Ψ i |Ψ i , ( 2.5)
where D = e r is the E1 operator, the subscripts i and f denote initial and final valence orbitals, respectively. Expanding the total wave function as given in Eq. (2.4) and retaining the terms only up to first order in G F , we get

E1 P NC = G F Ψ (0) f |D|Ψ (1) i + Ψ (1) f |D|Ψ (0) i Ψ (0) f |Ψ (0) f Ψ (0) i |Ψ (0) i = G F I =i Ψ (0) f |D|Ψ (0) I Ψ (0) I |H N SI AP V |Ψ (0) i E i -E I +G F J =f Ψ (0) f |H NSI AP V |Ψ (0) J Ψ (0) J |D|Ψ (0) i E f -E J , ( 2.6)
where the subscripts I and J represent the intermediate unperturbed states. We have used, here, the explicit form for the first order wave function given by

|Ψ (1) n = 1 G F I =n |Ψ (0) I Ψ (0) I |H NSI AP V |Ψ (0) n E n -E I . (2.7)
An important question we address in this paper is: How significant are the contributions from those states which were considered approximately in the sum-over-states approach and how they vary with the size of the systems? We would address this by carrying out a comparative study of E1 P NC results in two systems, namely Cs and Fr, of different atomic sizes. In order for the contributions of these higher excited states to be included, it is necessary to solve the first order perturbation equation directly. In other words, it is necessary to solve the equation following equation

(H DC -E (0) )|Ψ (1) n = 1 GF (E (1) -H NSI AP V )|Ψ (0) n , ( 2.8)
where E (1)

(= Ψ (0) n |H NSI AP V |Ψ (0)
n ) is the first order correction to E (0) which, however, vanishes in the present case.

III. APPLICATION OF RCC THEORY TO APV

The RCC method, which is equivalent to all order perturbation theory, has been used in the recent past and accurate results have been reported for many single valence systems [8,9,11,12]. In the RCC framework, the wave function of a single valence atom can be expressed as

|Ψ (0) n = = e T (0) {1 + S (0) n }|Φ n , ( 3.1)
where |Φ n is the reference state constructed from the Dirac-Fock (DF) wave function |Φ 0 of the closed-shell configuration by appending the valence electron n, that is, |Φ n = a † n |Φ 0 where a † n represents a creation operator which creates the valence electron n. Here T (0) and S (0) n are the RCC excitation operators which excite electrons from |Φ 0 and |Φ n , respectively, due to the residual Coulomb interactions. The corresponding excitation amplitudes are obtained by solving the following equations

Φ L |{ H DC N e T (0) }|Φ 0 = 0 (3.2) Φ L n |{ H DC N e T (0) }S n |Φ n = -Φ L n |{ H DC N e T (0) }|Φ n + Φ L n |S n |Φ n ∆E (0) n , ( 3.3)
with the superscript L(= 1, 2) representing the singly and doubly excited states from the corresponding reference states and the wide-hat symbol over H DC N e T (0) represent the linked terms of normal order atomic Hamiltonian H DC N and RCC operator T (0) . In the CCSD (CC with single and double excitations) approximation, the corresponding RCC operators are defined by

T (0) = T (0) 1 + T (0) 2 (3.4)
and

S (0) n = S (0) 1n + S (0) 2n (3.5)
for the closed-shell and single valence open-shell systems, respectively. The quantity ∆E (0) n in the above expression is the electron affinity energy (or negative of the ionization potential (IP)) for the valence electron which is evaluated by

∆E (0) n = Φ n |{ H DC N e T (0) }{1 + S (0) n }|Φ n . (3.6)
In addition to having considered full singles and doubles in the CCSD equations given in Eq. ( 3.3), we have also included the contributions from the important triple excitations perturbatively (known in the literature as CCSD(T) method) by defining

S (0),pert 3n = H DC N T (0) 2 + H DC N S (0) 2n , ( 3.7)
where the superscript pert denotes the perturbation and their contributions to ∆E

(0)
n are evaluating as

∆E (0),trip n = T (0) † 2 S (0),pert 3n . (3.8)
After solving for the amplitudes of T (0) , we solve Eqs. (3.3) and (3.6) simultaneously and obtain the amplitudes of S (0) n operators. Now, the total atomic wave function in the presence of H NSI AP V is expressed, in the RCC ansatz, as

|Ψ n = e T (0) +GF T (1) {1 + S (0) n + G F S (1) n }|Φ n , (3.9)
where T (1) and S

(1)

n are the first order perturbed amplitudes corresponding to the unperturbed RCC operators T (0) and S (0) n , respectively. On expanding the above equation keeping the terms only up to first order in G F yields

|Ψ n = e T (0) {1 + S (0) n + G F T (1) (1 + S (0) n ) + G F S (1) n }|Φ n . (3.10)
Comparing the above equation with Eq. (2.4), we get

|Ψ (1) n = e T (0)
{T (1) (1 + S (0) n ) + S (1) n }|Φ n . (3.11) In order to calculate |Ψ

(1) n as a solution of Eq. (2.8) in the RCC theory, we solve the excitation operator amplitudes of T (1) and S

(1) n using the following equations

Φ L |{ H DC N e T (0) T (1) }|Φ 0 = -Φ L | H NSI AP V e T (0) |Φ 0 (3.12)
and

Φ L n |{ H DC N e T (0) }S (1) n |Φ n = -Φ L n |{ H DC N e T (0) T (1) (1 + S (0) n ) + H NSI AP V e T (0) (1 + S (0) n )}|Φ n + Φ L n |S (1) n |Φ n ∆E n , ( 3.13)
after solving Eq. (3.2) and Eq. (3.3), respectively. In the above expression, notation H NSI AP V e T (0) is used for the connecting terms between H NSI AP V and T (0) . To keep the level of approximation uniform through out, both T (1) and S (1) n are truncated at single and double excitations by defining and S (1) 2n correspond to the perturbed double excitations, from closed-and open-shells, respectively. Since both the perturbed single and double excitation amplitudes are solved simultaneously, certain correlation effects due to the perturbed double excitations also reflect indirectly in the contributions of the perturbed single excitations.

T (1) = T (1) 1 + T (1) 2 (3.
After obtaining the unperturbed and the perturbed RCC operator amplitudes in both the closed-shell and onevalence open-shell atoms, we proceed to calculate the E1 P NC amplitude as

E1 P NC = G F Ψ (0) n |D|Ψ (1) n + Ψ (1) n |D|Ψ (0) n Ψ (0) n |Ψ (0) n = G F Φ n |{1 + S (0) † n }D{T (1)
1

(1 + S (0) n ) + S (1) 1n }|Φ n Φ n |{1 + S (0) † n }N 0 {1 + S (0) n }|Φ n +G F Φ n |{S (1) † n + (1 + S (0) † n )T (1) † }D{1 + S (0) n }|Φ n Φ n |{1 + S (0) † n }N 0 {1 + S (0) n }|Φ n , ( 3.16)
where we define D = (e T (0) † De T (0) ) and N 0 = e T (0) † e T (0) . The non-truncative series for D and N 0 are expanded using the Wick's generalized theorem and are truncated when the terms are below fifth order of the Coulomb interaction. The core-valence and valence correlation contributions are obtained from DT

(1) 1

and D(T

(1) 1 {1 + S (0) n } + S (1)
1n ), respectively, along with their conjugate terms.

Corrections due to the normalization of the wave functions are accounted by evaluating

Norm = Ψ (0) n |D|Ψ (1) n + Ψ (1) n |D|Ψ (0) n N n 1 + N n , (3.17)
where

N n = Φ n |{1 + S (0) † n }N 0 {1 + S (0)
n }|Φ n . Although, the CCSD(T) method described here, accounts for the contributions from the important unperturbed triple excitations it fails to include the direct triple excitation contributions to the E1 P NC calculations. To account for, at least, the lowest order direct triple excitation contributions (minimum up to fourth order in Coulomb interaction), we construct them with the open-shell RCC operators perturbatively as follows

S pqr,(0) nab = H DC N T (0) 2 + H DC N S (0) 2n n + a + b -p -q -r , ( 3.18)
and

S pqr,(1) nab = H DC N T (1) 2 + H DC N S (1) 2n n + a + b -p -q -r , ( 3.19)
where i is the single particle energy of an orbital i. These operators are finally considered as parts of S (0) n and S

(1) n in Eq. (III). This can be called as lo-CCSDvT approximation, which also accounts important lower order valence triple excitation effects in the final property calculations.

IV. RESULTS AND DISCUSSIONS

A. Orbitals generation

We have used Gaussian type functions

F GT O (r i ) = r nκ e -αir 2 i (4.1)
to construct the DF orbitals where α i is an arbitrary parameter which has to be chosen and r i represents a radial grid given by

r i = r 0 e h(i-1) -1 (4.2)
where the step size h is taken to be 0.03, the radial grid is increased up to i = 750, r 0 is 2 × 10 -6 in atomic units and n κ is the radial quantum number of the orbitals. Here α i s are chosen to satisfy the even tempering condition

α i = α 0 β i-1 (4.3)
and we have chosen different α 0 and β values for different symmetries (l) (known as even tempered basis) and they are given in Table I. The finite size of the nucleus in these systems is accounted by assuming a two-parameter Fermi-nuclear-charge distribution for evaluating the electron density over the nucleus as given by

ρ(r i ) = ρ 0 1 + e (ri-c)/a , ( 4.4)
where ρ 0 is the density for the point nuclei, c and a are the half-charge radius and skin thickness of the nucleus. These parameters are chosen as

a = 2.3/4(ln3) (4.5)
and

c = 5 3 r 2 rms - 7 3 a 2 π 2 , ( 4.6)
where r rms is the root mean square radius of the corresponding nuclei which is determined by r rms = 0.836A 1/3 + 0.570 (4.7)

in fm for the atomic mass A.

B. Results

In Table II, we give the E1 P NC amplitude results obtained from various calculations. From this work, we present the results using the DF, CCSD, CCSD(T) and lo-CCSDvT methods for both Cs and Fr. In the same table, we also compare our results with previously reported results using various many-body methods. Our Cs result matches reasonably well with the other calculations, but our Fr result differs significantly. We give below individual contributions from various RCC terms and express them in terms of correlation diagrams and level of excitations in order to facilitate the readers to understand the role of various correlation effects and intermediate states. Briefly, we discuss here the methods used in the other calculations. The difference between our DF and lo-CCSDvT results gives an idea about the amount of total correlation effects through the present method in these calculations. Considering the same DF wave functions, we have also employed CCSD and CCSD(T) methods where we find that the CCSD results are larger in magnitude than the CCSD(T) results. However, the E1 P NC result increases in magnitude for Cs in the lo-CCSDvT approximation, but it decreases for Fr; indicating that the triple excitation effects are stronger in Fr. From a comparison between these results, we observe that the dominant triple excitations effects arise through the CCSD(T) method.

C. Discussions

Here we discuss briefly different reported calculations of the above results at the DC approximation. About two decades ago, Blundell et al [13] had employed the linearized CCSD(T) method (LCCSD(T) method) to evaluate the unperturbed wave functions in Cs and then they had used a sum-over-states approach to evaluate the E1 P NC amplitude for the 6s 2 S 1/2 → 7s 2 S 1/2 transition in Cs. However, they had scaled their wave functions to fit the calculated energies of different states with the experimental results which reproduced many atomic properties quite accurately, but that does not show that the method they had used is capable of producing accurate ab initio results. Shabaev et al [15] have also obtained results using a CI method with a local form of the DF wave functions. Their Cs result matches with the result of Blundell et al. In our previous work on Cs [16], we had also calculated this quantity by considering the same α 0 and β parameters for orbitals of all symmetries (known as universal basis). With the new parameters, the Gaussian basis orbitals produce better wave functions in the nuclear region which are verified by studying the hyperfine interactions that will be reported elsewhere. The convergence of the RCC amplitudes are better than 10 -8 in the present case than 10 -6 in [16] and it gives a slightly different result. Dzuba et el have carried out a few calculations of these quantities using Brueckner orbitals using a Green function technique (Feynman diagram approach) that takes into account various classes of correlation effects to all orders and avoids the sum-over-states approach [14,17,19]. Their results also differ from each other and in some cases with others as can be seen in Table I. The most recent calculation on Cs is reported by Porsev et al [18] using the RCC method that includes all single and double excitations with all valence triple excitations (CCSDvT method). However, they have finally used a sum-over-states approach to calculate the E1 P NC amplitude of the 6s 2 S 1/2 → 7s 2 S 1/2 transition. They give contributions from 6p 1/2 to 9p 1/2 singly excited states as "Main", which contributes -0.8823 (18), and the remaining contributions as "Tail", which is obtained using other many-body methods as -0.0175 (18), (results are always given in ×10 -11 iea 0 (-Q W /N ) here onwards) at the DC approximation. In contrast to their approach, we have used the lo-CCSDvT method, but have included contributions from the core and doubly excited states in a manner similar to that of the singly excited states. In Table III, we present the individual contributions from different RCC terms in our CCSD(T) method to the E1 P NC calculations for the 6s 2 S 1/2 → 7s 2 S 1/2 and 7s 2 S 1/2 → 8s 2 S 1/2 transitions in Cs and Fr, respectively. It is evident that the trends of the contributions from the different RCC terms are similar for both Cs and Fr. The important contributions come from three terms:

1/2 → 7s 2 S (0) 1/2 ) ( 6 s 2 S (0) 1/2 → 7s 2 S (1) 1/2 ) DF result DH NSI AP V -1.0460 H NSI AP V D 0.3224 T (1) † 1 D 0.0374 DT (1) -0.0355 D1S (1) 1i -1.7467 S (1) † 1f D1 0.1841 D1S (1) 2i -0.0046 S (1) † 2f D1 0.0226 S (0) † 1f D1S (1) 1i 0.4061 S (1) † 1f D1S (0) 1i 0.2141 S (0) † 2f D1S (1) 1i -0.0055 S (1) † 1f D1S (0) 2i 0.0301 S (0) † 1f D1S (1) 2i -0.0026 S (1) † 2f D1S (0) 1i -0.0028 S (0) † 2f D1S (1) 2i 0.0011 S (1) † 2f D1S (0) 2i -0.
1/2 → 8s 2 S (0) 1/2 ) ( 7 s 2 S (0) 1/2 → 8s 2 S (1) 1/2 ) DF result DH NSI AP V -19.4661 H NSI AP V D 5.9836 T (1) † 1 D 0.7147 DT (1) -0.6876 D1S (1) 1i -31.7107 S (1) † 1f D1 4.0538 D1S (1) 2i -0.0314 S (1) † 2f D1 0.5930 S (0) † 1f D1S (1) 1i 7.8001 S (1) † 1f D1S (0) 1i 4.1072 S (0) † 2f D1S (1) 1i -0.0423 S (1) † 1f D1S (0) 2i 0.6580 S (0) † 1f D1S (1) 2i -0.0615 S (1) † 2f D1S (0) 1i -0.0911 S (0) † 2f D1S (1)
D 1 S (1) 1i , D 1 S (1)
2i and DT and their corresponding conjugate terms where D 1 is the effective one-body terms of D and bare operator D is its lowest order term. Other terms correspond to higher order RCC terms, but they are not small. Contributions given as Others represent the RCC terms that come from the effective two-body terms of D after contracting with the open-shell RCC terms. We give diagrammatic representations of the above three RCC terms in Fig 2 (without their conjugate terms) and their lowest order terms. From this figure, it can be noticed that DT 1i contain the lowest order DF contributions from the core (hole) and virtual orbitals, respectively; hence D 1 S (1) 1i always gives the largest contribution. Again, the perturbed states arising through the ground state contribute predominantly while contributions from the perturbed excited states are comparatively smaller than those corresponding to the ground state but with the opposite signs. The final results are the outcome of these cancellations. The most core correlation effects are coming from DT (1) 1 and its conjugate terms while there are small contributions that come from DT and its conjugate terms which are included in "Others". As seen these correlation effects mostly cancel out with their corresponding conjugate terms in both the systems. The contributions from singly excited states can be estimated by summing the total contributions from

D 1 S (1) 1i , S (0) † 1f D 1 S (1) 1i , S (0) † 2f D 1 S (1)
1i and their conjugate RCC terms. The remaining contributions are come from the doubly excited states apart from the normalization corrections to the RCC wave functions. After accounting the corresponding corrections from the normalization of the wave functions to various RCC operators, we get contributions as 0.0018, -0.8980 and 0.0073 from the core correlation, singly excited states and doubly excited states, respectively, in Cs. Similarly, we get the contributions as 0.0271, -14.7551 and 0.3158 from the core correlation, singly excited states and doubly excited states, respectively, in Fr. Clearly, the core correlation and doubly excited states contributions in Fr are large. In Cs and Fr, the doubly excited states contributions are around 1% and 2%, respectively, with opposite signs. Hence, these contributions should also be accounted accurately for high precision results.

To understand more clearly about the role of different excited states that play important roles in the E1 P NC result of Cs through our RCC approach which have also been studied by other methods, we consider some of the important RCC terms and analyze their contributions through various excitations level. These results are given explicitly in Table IV. As noticed, only 4p 1/2 and 5p 1/2 core orbitals are important to account the core correlation effects; this is because the single particle energy levels between these orbitals with the valence electrons are small. The singly excited states from 6p 1/2 to 9p 1/2 contribute predominantly, but there are also significant contributions from some of

D T P S 1f (1) 1i S 2i (1)
(1)

D D D D P i i i f f f f f f f f i i i i i V V V V V
V i i i i i i i i i f f f f f f f f f V V V V V V V T 2 (0) S (1) S (1) T 2 (0) D T 1 (0) (1) 1

S

S

D T 1 (0)
T 1

T 1 FIG. 3: Some of the important contributing diagrams from the triple excitations which contribute significantly to the lo-CCSDvT method along with their conjugate terms.

the continuum states. The continuum orbital, 10p 1/2 , whose density in the nuclear region is large. There are various doubly excited states that play important roles, mainly from the 6skp to 7sld and 7skp to 6sld excitations for any arbitrary principal quantum numbers k and l, in calculating the above quantities but mostly they cancel each other. Therefore, a suitable choice of basis functions and an accurate many-body method are necessary for accounting these contributions precisely.

To compare the role of the correlation effects in calculating E1 P NC amplitudes in Cs and Fr, we compare the DF results and then their final results. It is about 20% in Cs, but it is only about 7% in Fr. The small correlation effects in Fr is because of strong cancellations between the contributions from the initial and final perturbed states. Again, it is interesting to note that the CCSD results were larger than the CCSD(T) results suggesting that there were cancellations from the triple excitations. Using the lo-CCSDvT method, these results increase in Cs slightly but cancel out in Fr. The triple excitation effects in Fr are strong enough to be taken seriously into account for high precision results. In Fig. 3, we give some of the important RCC diagrams (without their conjugate terms) that give large contributions through the lo-CCSDvT method. In the future, we would like to study other relevant properties in order to determine the accuracy of our calculated E1 P NC amplitudes in these systems.

V. CONCLUSION

We have applied the relativistic coupled-cluster method to calculate the atomic parity violating effects in Cs and Fr which involve the interplay of the long range electrostatic interactions and the short range weak interactions. The present approach considers the Coulomb interaction up to all orders and the weak interaction between the nucleus and the electrons up to first order. From the detailed analysis of different contributions, it is clear that the doubly excited states play an important role in obtaining precise results. These contributions are quite large in magnitude, however, with opposite signs, than those of the singly excited states in Fr, resulting in a large cancellation in the final results. They are about 1% in Cs, for which the experimental accuracies have been claimed to be around 0.35%. This suggests that, it requires a method like our present approach to consider them accurately for high precision calculations.

FIG. 1 :

 1 FIG. 1: Diagrammatic representation of the electron-nucleus interactions due to the electromagnetic (with the exchange of photon γ denoted by solid curved line) and the weak (with the exchange of heavy boson Z0 denoted by straight dashed line) interactions.

  correspond to the perturbed single excitations and T (1) 2

FIG. 2 :

 2 FIG.2: Break-down of the perturbed RCC diagrams into the lower-order many-body theory (MBPT) diagrams. Here i and f represent the initial and final orbitals, P is the PNC operator and V is the two-body part of the Coulomb interaction Hamiltonian. Lines going up and down with single arrows represent virtual and hole (core) orbitals and lines with double arrows represent attached valence orbitals to the closed-shell DF wave function.

TABLE I :

 I Used α0 and β parameters for different symmetries (l) to construct DF orbitals using GTOs in Cs and Fr.

	H H H H H α0/β l	0	1	2	3	4
	α0	0.00190 0.001825 0.00183 0.00185 0.00187
	β	2.91	2.94	2.94	3.05	3.09

TABLE II :

 II Comparison of E1P NC results of 133 Cs and 223 Fr due to electron-electron Coulomb interactions from various methods in ×10 -11 iea0(-QW /N ). S 1/2 → 7s 2 S 1/2 7s 2 S 1/2 → 8s 2 S 1/2

	133 Cs	223 Fr	Method	Reference
	6s 2 -0.7236	-13.4825	DF	This work
	-0.8941	-14.5134	CCSD	This work
	-0.8889	-14.4219	CCSD(T)	This work
	-0.8892	-14.4106	lo-CCSDvT	This work
	-0.904(9)		LCCSD(T) + SL + SS	[13]
	-0.9078		BO + GFCP	[14]
	-0.904	-15.72	CI	[15]
	-0.902(4)		CCSD(T)	[16]
	-0.9001	-15.229	BO + GFCP	[17]
	-0.8998(25)		CCSDvT + SS	[18]
		-15.9	BO + GFCP	[19]
		-15.41(17)	LCCSD + RPA + ExpEng + SS	[20]
	NOTE: Contributions from Breit interaction, QED correction and nuclear effects are not considered here.
	Abbreviations SS	: sum-over-states
	SL	: scaling
	BO	: Brueckner orbitals
	GFCP : Green function technique for all order correlation potential
	CI	: configuration interaction method
	ExpEng : experimental energy
	RPA	: random phase approximation

TABLE III :

 III Contributions to the E1P NC calculations in Cs and Fr using CCSD(T) method (in ×10 -11 iea0(-QW /N )).

	(A) Cs
	Initial pert. terms	Final pert. terms
	(6s 2 S	(1)

TABLE IV :

 IV Important contributions to the E1PNC calculations in Cs from various intermediate states (I and J notations are used as per Eq. (2.6) in ×10 -11 iea0(-QW /N )).

	6s 2 S	(1) 1/2 → 7s 2 S	(0) 1/2	6s 2 S	(0) 1/2 → 7s 2 S	(1) 1/2
							I	Results	J	Results
	T	(1) † 1	D	D T (1)
							4p 1/2 → 7s	0.0005	4p 1/2 → 6s	-0.0005
							5p 1/2 → 7s	0.0368	5p 1/2 → 6s	-0.0349
	D1S	(1) 1i	S	(1) † 1f D1
							6s → 6p 1/2	-1.9195	7s → 6p 1/2	1.7648
							6s → 7p 1/2	0.1532	7s → 7p 1/2	-1.5022
							6s → 8p 1/2	0.0347	7s → 8p 1/2	-0.0902
							6s → 9p 1/2	-0.0084	7s → 10p 1/2	0.0093
							6s → 10p 1/2	-0.0051
	S	(0) † 1f D1S	(1) 1i	S	(1) † 1f D1S	(0) 1i
							6s → 6p 1/2	0.1716	7s → 6p 1/2	0.2556
							6s → 7p 1/2	0.1302	7s → 7p 1/2	0.1213
							6s → 8p 1/2	0.0258	7s → 8p 1/2	-0.0611
							6s → 9p 1/2	0.0735	7s → 9p 1/2	-0.0946
							6s → 10p 1/2	0.0053	7s → 10p 1/2	-0.0055
	D1S	(1) 2i

VI. ACKNOWLEDGMENT

I am grateful to B. P. Das for useful discussions and H. S. Nataraj for his critical reading of the manuscript. This work was supported by NWO under the VENI program with Project No. 680-47-128. I thank the C-DAC TeraFlop Super Computing facility, Bangalore, India for the cooperation to carry out these calculations on its computers.

S (1) † 2f D1 6s5p 1/2 → 7s6s -0.0019 7s5p 3/2 → 6s9s 0.0002 6s5p 1/2 → 7s8s -0.0001 7s5p 3/2 → 6s10s 0.0009 6s5p 3/2 → 7s8s 0.0001 7s5p 1/2 → 6s6d 3/2 0.0001 6s5p 1/2 → 7s9s -0.0009 7s5p 3/2 → 6s6d 3/2 -0.0003 6s5p 3/2 → 7s9s 0.0009 7s5p 3/2 → 6s7d 3/2 -0.0001 6s5p 1/2 → 7s10s -0.0014 7s5p 1/2 → 6s8d 3/2 0.0006 6s5p 3/2 → 7s10s 0.0001 7s5p 3/2 → 6s8d 3/2 -0.0003 6s5s → 7s9p 1/2 0.0001 7s5p 1/2 → 6s9d 3/2 0.0019 6s5s → 7s10p 1/2 0.0002 7s5p 3/2 → 6s9d 3/2 0.0003 6s5p 1/2 → 7s6d 3/2 0.0006 7s5p 1/2 → 6s10d 3/2 0.0005 6s5p 3/2 → 7s6d 3/2 0.0005 7s5p 3/2 → 6s10d 3/2 0.0001 6s5p 1/2 → 7s7d 3/2 0.0003 7s5p 3/2 → 6s5d 5/2 0.0044 6s5p 3/2 → 7s7d 3/2 0.0002 7s5p 3/2 → 6s6d 5/2 0.0018 6s5p 1/2 → 7s8d 3/2 0.0013 7s5p 3/2 → 6s7d 5/2 0.0008 6s5p 3/2 → 7s8d 3/2 0.0011 7s5p 3/2 → 6s8d 5/2 0.0041 6s5p 1/2 → 7s9d 3/2 0.0006 7s5p 3/2 → 6s9d 5/2 0.0048 6s5p 3/2 → 7s9d 3/2 0.0005 7s5p 3/2 → 6s10d 5/2 0.0003 6s5p 1/2 → 7s10d 3/2 -0.0004 7s4d 3/2 → 6s8f 5/2 0.0003 6s5p 3/2 → 7s10d 3/2 -0.0002 7s4d 3/2 → 6s9f 5/2 0.0005 6s5p 3/2 → 7s8d 5/2 -0.0003 7s4d 5/2 → 6s8f 7/2 0.0004 6s5p 3/2 → 7s9d 5/2 -0.0013 7s4d 5/2 → 6s8f 7/2 0.0004 6s5p 3/2 → 7s10d 5/2 -0.0012 7s4d 5/2 → 6s9f 7/2 0.0009 6s4d 3/2 → 7s8f 5/2 -0.0001 6s4d 3/2 → 7s9f 5/2 -0.0004 6s4d 5/2 → 7s8f 7/2 -0.0005