J R Armstrong

T Zinner

D V Fedorov

A S Jensen

Analytic Harmonic Approach to the N -body problem

Keywords: numbers: 03.75.Hh, 05.30.Jp, 21.45.+v

We consider an analytic way to make the interacting N -body problem tractable by using harmonic oscillators in place of the relevant two-body interactions. The two-body terms of the N -body Hamiltonian are approximated by considering the energy spectrum and radius of the relevant two-body problem which gives frequency, center position, and zero point energy of the corresponding harmonic oscillator. Adding external harmonic one-body terms, we proceed to solve the full quantum mechanical N -body problem analytically for arbitrary masses. Energy eigenvalues, eigenmodes, and correlation functions like density matrices can then be computed analytically. As a first application of our formalism, we consider the N -boson problem in two-and three dimensions where we fit the two-body interactions to agree with the well-known zerorange model for two particles in a harmonic trap. Subsequently, condensate fractions, spectra, radii, and eigenmodes are discussed as function of dimension, boson number N , and scattering length obtained in the zero-range model. We find that energies, radii, and condensate fraction increase with scattering length as well as boson number, while radii decrease with increasing boson number. Our formalism is completely general and can also be applied to fermions, Bose-Fermi mixtures, and to more exotic geometries.

Introduction

Determining the ground state and properties of N interacting particles in some fixed geometry is at the core of many disciplines in physics and other natural sciences. However, in general even for moderate values of N , methods based on first principles are either intractable or extremely time-consuming. Fortunately, the properties of many systems can be described by correlations that involve just a few particles and the problem of many particles can be reduced to consider a much smaller number in the background of the remaining particles. Two-body interactions still have a tendency to make even such few-body problems very difficult to solve and insights gained from approximations that allow analytic treatments are therefore very useful.

In configuration interaction methods, where large basis states of properly symmetrized wave functions are built and diagonalized to determine system properties, a basis of harmonic oscillator states can be very convenient as matrix elements of the two-body interaction are easy to calculate. This approach was used early on in the context of the nuclear shell model [START_REF] Jensen | Elementary Theory of Nuclear Shell Structure[END_REF][START_REF] Heyde | The Nuclear Shell Model[END_REF]. Turning this method upside down the potentially complicated two-body interaction could be reproduced or simulated by a simple harmonic oscillator potential. The great advantage is obviously that the coupled set of differential equations of motion are much simpler to solve, and a number of properties are easily systematically obtained as functions of interaction parameters and particle number.

In subfields of physics where the structure is already established at a given level of accuracy, the insight gained from oscillator approximations is most often insufficient for improvements. However, in cold atomic gas physics it is now possible to prepare systems with very exotic and often unknown properties, to use controlled two-body interactions, to study different geometries and dimensions, and to vary trapping conditions [START_REF] Bloch | [END_REF]. Analytical oscillator approximations can therefore be expected to be very valuable, as it has been in other fields of physics. While we are concerned mainly with the quantum mechanical many-body problem here, we note that similar techniques are currently used also in the study of classical dynamics [4].

Obviously, the accuracy of the oscillator approximation increases as the potentials resemble oscillators. This means that pronounced smooth potential single minima with room for bound states are directly suited for investigations of system structures as functions of particle number and other characteristics. However, also much weaker binding potentials would reveal correct overall qualitative, and perhaps also semiquantitative, behaviour in an appropriate oscillator approximation. This is especially emphasized by the universal behaviour of a number of weakly bound structures which only depend on integral properties like scattering length. In that case, the bulk part of the potentials are not crucial by themselves but rather the large distance effect or equivalently the tail of the wave function or the binding energy. All these quantities are related in the correct model descriptions but for specific purposes a subset may suffice. It seems clear that continuum properties like scattering behaviour are beyond the regime where useful results can be expected. Still even phase shifts can be extracted from properly discretized continuum states [5].

The harmonic oscillator has been used in virtually every aspect of physics where potentials are needed. Still, in full generality, the procedure is not well described probably for several reasons. First, relative and center-of-mass (cm) motions are not separable even for two interacting particles with different masses in a trapping one-body potential. Second, the simplest approximation for a self bound N -body system is found by the mean-field approximation where the spurious cm-motion is ignored. Third, oneor two-body potentials centered at different points in space have only been of interest when the harmonic approximation is insufficient, as e.g., in chemistry and for crystal structures. Now optical lattices and split traps offer the possibility to also use multicentered potentials. If only the relative motion is of interest, a decoupling scheme is necessary in order to separate the full solution to relative and cm-motions.

For cold atoms the oscillator approximation has been applied recently in [6,7,8]. These works considered Bose-condensate properties for equal mass particles, and this is almost the only case where the center-of-mass motion separates. In this report we develop formalism to treat the most general harmonically interacting system for bosons (or distinguishable particles). In turn, we solve the N -body problem for arbitrary quadratic forms of the one-and two-body interactions. The use of Cartesian coordinates allows simple solutions for both one, two and three spatial dimensions, and any nonspherical behaviour of the corresponding potentials. First we derive the transformation matrices from initial to final particle coordinates where the differential equations are completely decoupled. We also derive expressions for various quantities like energies, root-mean-square radii, density matrix and its dominating eigenvalue. The hope is to provide simple tools to reveal at least qualitative features of the new and unknown systems under design and investigations particularly in the field of cold atomic gases.

We demonstrate our method by finding solutions in two (2D) and three (3D) dimensions for N identical and pairwise interacting bosons in external harmonic onebody potentials. The pairwise interactions are taken from the celebrated results of Busch et al. [9] for two particles with zero-range interaction in a harmonic trap, the validity of which have been tested in ultracold atomic gas experiments [10]. Our method thus provides an analytical approximation to the N -boson problem with short-range interactions. The condensate fraction is readily available from our calculations and shows interesting behaviour as the scattering length is tuned and the number of particles changes. In particular, at small positive scattering length we find a highly fragmented state in both two and three dimensions.

Theoretical derivations

We first define the Hamiltonian in its most general quadratic form for both coordinates and kinetic energy derivative operators. The spin degrees of freedom are omitted and symmetries of the spatial wave functions can in principle easily be imposed by permutations of the coordinates. We use matrices to simplify the derivations. We then derive the coordinate transformation to decouple the set of coupled oscillators and distinguish between cases where the cm motion is free and confined by one-body potentials. Lastly, we calculate pertinent properties.

Hamiltonian

We consider a system of N , possibly different, particles of mass m k (k = 1, . . . , N) interacting through deformed harmonic potentials V int . The particles are in addition subject to external one-body potentials, V ext , for each particle constructed as a sum of harmonic oscillators with different centers. The total Hamiltonian H = T + V with kinetic energy T and potential V = V int + V ext is then given by

T = - N k=1 2 2m k ∂ 2 ∂x 2 k + ∂ 2 ∂y 2 k + ∂ 2 ∂z 2 k , (1)
V int = 1 4 N i,k=1 V ik,0 + µ ik ω 2 x,ik (x i -x k + x ik,0 ) 2 (2) + ω 2 y,ik (y i -y k + y ik,0 ) 2 + ω 2 z,ik (z i -z k + z ik,0 ) 2 , V ext = 1 2 N k=1 m k ω 2 x,k (x k -x k,0 ) 2 (3)
+ ω 2 y,k (y k -y k,0 ) 2 + ω 2 z,k (z k -z k,0 ) 2
, where (x k , y k , z k ) are the (x, y, z)-coordinates of the k th particle, µ ik = m i m k /(m i +m k ) is the reduced mass of particles i and k, (ω x,ik , ω y,ik , ω z,ik ) are the frequencies in the (x, y, z)-directions for the interaction potential between particles i and k, and (ω x,k , ω y,k , ω z,k ) are the frequencies in the (x, y, z)-directions for the one-body potential on the k th particle with centers specified by (x k,0 , y k,0 , z k,0 ). The factor 1/4 is made of two factors 1/2 where one of them is the conventional notation for an oscillator potential. The other factor 1/2 is to count the two-body interaction only once when the i, k summations are extended to assume all integer values from 1 to N . The shift of the interaction centers, x ik,0 , for each pair of the two-body interactions should then change sign when the particles are interchanged, x ik,0 = -x ki,0 , which implies that the diagonal has to vanish, x ii,0 = 0, in accordance with zero self interaction.

This Hamiltonian has the most general quadratic form expressed in terms of the parameters for one-and two-body oscillator potentials. The shift of potential energy, V ik,0 , of the energy for each pair allows adjustment without change of structure. The shifts of both one-and two-body oscillator centers suggest applications approximating optical lattice potentials. The frequencies traditionally all enter as squares which suggest attraction but the formalism is equally applicable for imaginary frequencies or equivalently negative values of these squared frequencies. To produce stable solutions with such repulsive interactions requires sufficient attraction from the other two-body interactions or from the external fields. The choice of Cartesian coordinates allows independent solution for each dimension, and thereby treats deformations and different dimensions without any additional complications. Obviously this also prohibits direct use of symmetries and conserved quantum numbers where the dimensions are mixed. One example is angular momentum conservation in the absence of external fields.

We now proceed by rewriting the Hamiltonian in matrix form. For this we define vectors, x = (x 1 , x 2 , ..., x N ) T , ∇ x = (∂/∂x 1 , ∂/∂x 2 , ..., ∂/∂x N ) T , where a vector is given as a column of its coordinates, which means the transposed of the row as indicated by "T ". The y and z-direction can be defined analogously if necessary. We treat each coordinate independently and may therefore omit the x-index to simplify the notation in the derivation. The x-part of the Hamiltonian, H x , in (1)-( 3) is given by:

H x = 1 2 ∇ T x T ∇ x + 1 2 x T A x + c • x + V shift , (4)
where the kinetic energy matrix, T ik = -δ ik 2 /(m i ), is diagonal and depends only on inverse masses. The constant term, V shift consists of the sum of all separate shift energies:

V shift = 1 2 N k=1 m k ω 2 x,k x 2 k,0 (5)
+ 1 4 N i,k=1 V ik,0 + µ ik ω 2 x,ik x 2 ik,0
.

The quadratic potential term in (4) contains the symmetric matrix A which is given in terms of masses and frequencies by

A i =k = -µ ik ω 2 x,ik (6)
A kk = m k ω 2 x,k + N i,i =k µ ik ω 2 x,ik . (7)
The components of the coefficient vector c in the linear term are

c k = -m k ω 2 k x k,0 + N i=1 µ ik ω 2 x,ik x ik,0 . (8)
The y and z-parts of the Hamiltonian, H, are completely analogous and we have

H = H x + H y + H z .

Reduction to standard form

The linear term in (4) can be eliminated by translating the coordinates by

x = x -a , x = x + a , (9)
where the translation vector a is determined by the requirement that all terms linear in x k must vanish from the Hamiltonian. This condition amounts to A a = -c. As containing only second derivatives the kinetic energy operator remains unchanged by this linear translation. In total the Hamiltonian in the new coordinates has only the same quadratic terms in both coordinates and derivatives. However, linear terms have disappeared and the term, -1/2 a T A a, should be added to V shift in (5). Any solution, a, to A a = -c eliminates the linear terms. Many solutions always can be found, but a unique solution only exists when A is non-singular, that is A is invertible. A singular A-matrix is equivalent to a subset of non-interacting particles which also all are unaffected by the one-body potentials. Then the corresponding degrees of freedom are already decoupled and following the trivial motion in free space. As already decoupled they should therefore from the beginning be removed from the reduction procedure. The dimension of the A-matrix is correspondingly reduced.

In practice, only one frequently occurring example is interesting and requiring special attention. This is when all interactions are translation invariant and only depending on coordinate differences between particles. Then the center-of-mass motion is as in free space. External one-body fields acting on all particles act on the center-ofmass coordinate and remove the corresponding degeneracy.

Previous works with oscillators [6,7,8] considered equal mass systems where separation of the center-of-mass motion is straightforward. Here we provide a general procedure valid for all sets of masses and all one-and two-body interactions. In all cases we transform to relative and center-of-mass, X, coordinates, where the latter is defined by

XM = N k=1 m k x k , M = N k=1 m k , ( 10
)
where M is the total mass of the system. Choosing the new set of coordinates, x, by xi ≡ x i -X, for i = 1, 2, ..., N -1, supplemented by xN ≡ X, we define the transformation matrix, F , by

x T = F xT , (11)
where the transformation matrix takes the form

F =      1 0 . . . 1 0 1 . . . 1 . . . . . . . . . 1 -m 1 m N -m 2 m N . . . 1      . ( 12
)
The specific labelling singles out the N th (last) particle with mass, m N , in the transformation F . The inverse transformation, F -1 , from original to new coordinates is explicitly given by

F -1 =      1 -m 1 M -m 2 M . . . -m N M -m 1 M 1 -m 2 M . . . -m N M . . . . . . . . . -m N M + m 1 M + m 2 M . . . + m N M      . ( 13
)
The Hamiltonian is now easily transformed to the new set of coordinates, by direct insertion of (12) for the coordinates and the inverse and transposed, (F -1 ) T , for the derivative kinetic energy operators, ∇ x. The kinetic energy immediately separates into a sum of relative and center-of-mass dependent terms, whereas the potential part separates when it is translationally invariant, and otherwise not. In any case these coordinates can be used in the following derivations.

Decoupling the oscillators

We assume that we have the quadratic form without linear terms. We transform to relative and center-of-mass coordinates in all cases even though we could have worked with the initial coordinates. We rename the coordinates and omit the primes in the following expressions. We begin by diagonalizing the quadratic part of the potential, which then is the F T AF matrix (see ( 4), ( 6) and ( 7)). The orthonormal coordinate transformation Q corresponding to the diagonalization is defined by the requirements

Q T F T AF Q = D , x = F x = F Q t , ( 14
)
where D is diagonal with eigenvalues d i .

If A is singular at least one of the eigenvalues, d i , is zero, as when the interactions only depend on relative coordinates. One of these zero eigenvalues then has an eigenvector corresponding to the center-of-mass coordinate. This mode now has to be fully decoupled after the kinetic energy is expressed in the same coordinates. This is achieved in general by replacing the zero eigenvalue by any finite number. If only relative motion is of interest the Hilbert space spanned by this eigenvector should simply be removed.

We proceed to perform a new non-orthonormal transformation of the coordinates defined by t = D u, where D is diagonal and given by Dik = δ ik d 0 /d i . The number d 0 can be chosen arbitrarily. We choose to maintain the total norm which implies that Π N i=1 d

(x) i = d N 0 .
This transformation is nothing but scaling the lengths of each of the eigenvectors to let D become proportional to the unit matrix. The transformation from x to u-space of the derivative vector is then given as ∇ x = (F -1 ) T Q( D-1 ) T ∇ u . The kinetic energy matrix in the u-coordinates is D-1 Q T (F -1 )T (F -1 ) T Q( D-1 ) T , where the space corresponding to the center-of-mass coordinate decouples from all the other degrees of freedom. This kinetic energy matrix is diagonalized, i.e.

T x = 1 2 ∇ T u D-1 Q T F -1 T (F -1 ) T Q( D-1 ) T ∇ u = 1 2 ∇ T v P T D-1 Q T F -1 T (F -1 ) T Q( D-1 ) T P ∇ v = 1 2 ∇ T v T ∇ v , (15)
where the orthonormal transformation

∇ u = P ∇ v , or u = P T v, is chosen to make T diagonal, i.e. Tkk = -2 /( mk ), or equivalently diagonalize D-1 Q T F -1 T (F -1 ) T Q( D-1 ) T .
This final orthonormal transformation leaves the potential energy part as a diagonal matrix with unchanged diagonal elements, i.e.

V x = 1 2 x T F T AF x = 1 2 t T D t = 1 2 u T DT D D u = 1 2 d 0 u T u = 1 2 d 0 v T v = 1 2 v T Vx v , ( 16
)
where Vx is the diagonal unit matrix with elements d 0 = mk ω2 k , which defines the output frequencies, ωk . This is a factorization return to the ordinary oscillator notation by using the masses determined from the kinetic energy eigenvalues in (15). If in the process of transforming from the original xcordinates to the v coordinatees zero eigenvalues were generated, then care must be taken to make sure sure that they do not contribute to any final observable. This is easily achieved with ωk = 0 for these modes.

The total Hamiltonian for the x-coordinate is now a set of decoupled oscillators in the new coordinates v, i.e.

H x = V shift (17)
+ N k=1 - 2 2 mx,k ∂ 2 ∂v 2 x,k + 1 2 mx,k ω2 x,k v 2 x,k ,
where we inserted the label x on masses, oscillator parameter and v-coordinates. The harmonic oscillator eigenvalues and eigenfunctions correspond to the frequencies ωx,k and the length parameter b x,k as usual given by b 2

x,k = /( mx,k ωx,k ). In total the transformation M from initial to new coordinates and vice versa are

v = M( x -a) , x -a = M -1 v , (18)
M = D-1 Q T F -1 . ( 19
)
The normal mode of the system is expressed by the eigenvector, and the corresponding eigenvalue indicates the ease or difficulty of exciting that particular normal mode.

Basic properties

Observables expressed as expectation values of operators O are found by

Ψ|O|Ψ = d N xd N yd N z (20) Ψ * ( v x , v y , v z )O( x, y, z)Ψ( v x , v y , v z ) ,
where the wave functions are simplest in the transformed coordinates while the operators probably are simpler in the original particle coordinates. We have only considered the spatial wave function Ψ. Effects of spin dependence of interaction or quantum statistics have to be inserted separately.

The wave functions are products of the three one-dimensional harmonic oscillator wave functions in the new coordinates, that is Gaussians, exp(-v 2

x,k /(2b 2 x,k )), times Hermite polynomials with arguments v x,k /b x,k for each contributing mode k. Other analogous products arise from the y and z-directions. The normalization is as usual when we use the final (v x,k ) coordinates in the wave function. All the above transformations, including that of the non-orthonormal matrix, D, was chosen as total norm conserving. Therefore the volume elements have the same structure in initial and new coordinates, i.e. Π N k=1 dx k = Π N k=1 dv x,k . The expectation value of the Hamiltonian gives the energy, which for the eigenmode, k, of given quantum number, n x,k , is ωx,k (n x,k + 1/2). Adding the similarly obtained results from the y and z-direction we get the total energy for a given set, (n x,k , n y,k , n z,k ), of quantum numbers

E n x,k ,n y,k ,n z,k = V shift + N k=1 (ω x,k (n x,k + 1/2) (21) + ωy,k (n y,k + 1/2) + ωz,k (n z,k + 1/2)) ,
where the frequencies of the non-contributing modes are inserted as zero.

Beside energies the simplest observables are sizes, i.e. average values of the interparticle distance(s) in the system. The spatial extension of the probability for each particle is measured by its root mean square radius. With the relation in ( 18) we find for the ground state Ψ

Ψ|x 2 i |Ψ = Ψ| (A -1 c) i + N k=1 (F Q DP T ) ik v k 2 |Ψ = ((A -1 c) i ) 2 + N k=1 1 2 b 2 x,k ((F Q DP T ) ik ) 2 , (22)
where we used that the odd powers of the coordinates vanish after integration. If the center of mass coordinate X is decoupled we exclude that term in the summation, and the resulting sum is the computation of the expectation of (x i -X) 2 . The other dimensions should be added.

One-body density matrix

The simplest correlation function is the one-body density matrix, ρ(x 1 , x 1 ), which is the starting point of the calculation of statistical properties. It is interesting in itself but serves here also as a simple illustration of analytical calculations. The one-body density matrix is relevant for bosons (or distinguishable particles) which is our main application of the method to be discussed in the subsequent section. For fermions, the one-body density matrix is a rather uninteresting object due to the Pauli principle and off-diagonal long-range order (as arising for instance from pairing) is determined by the two-body density matrix [11]. While we do not consider fermions explicitly in the current presentation, we address the general issue of symmetrization in the subsection below for completeness. Note here that we have taken great care to remove the center-of-mass coordinate so that only intrinsic coherence is studied. This is particularly important if N is very small [12] or if the system is subjected to a temporally or spatially varying additional potential [13].

The ground state wave function is a product of Gaussians in the new coordinates, v with different units of length b k . The exponent can be written -v T B v where the matrix B is diagonal with elements B kk = 1/(2b 2 k ). Using (18) we return to the initial coordinate where

ln Ψ = -v T B v = x T Z x + constant , (23)
which defines the matrix Z when we assume all interaction centers are at the origin. The exponent of the density matrix for identical bosons, where we select the particle labeled 1, is then

-x 1 Z 11 x 1 -x 1 Z 11 x 1 (24)
-

(x 1 + x 1 ) N k=2 (Z 1k + Z k1 )x k + 2 N i,k=2 x i Z ik x k .
To integrate over all other coordinates (from -∞ to +∞) than x 1 and x 1 we complete the squares. We define the vector w = Z 1k + Z k1 , for 2 ≤ k ≤ N and make the substitution x = f -q, where q = 1/2(x 1 + x 1 )( Z + ZT ) -1 w, and Z is the Z matrix without the first row and column. After integration we are left with the density matrix

ρ(x 1 , x 1 ) = N exp [-(x 1 Z 11 x 1 + x 1 Z 11 x 1 ) ( 2 5 )
+ 1 4 (x 1 + x 1 ) 2 w T ( Z + ZT ) -1 w] ,
where the normalization, N is

N = π(Z 11 - 1 2 w T ( Z + ZT ) -1 w). (26)
Going further, we can re-write the exponent as

- 1 2 Z 11 (x 1 -x 1 ) 2 (27)
+ 1 4 w T ( Z + ZT ) -1 w - 1 2 Z 11 (x 1 + x 1 ) 2 ,
where the ratio d x is given as

d x = 2Z 11 2Z 11 -w T ( Z + ZT ) -1 w . ( 28
)
This ratio determines the largest eigenvalue, λ, obtained after diagonalization of the density matrix, i.e.

λ = 2 1 + √ d x , (29)
where the subscript x is to remind us that this expression is valid for one dimension.

In more dimensions the wave functions are products and consequently also the density matrix and its eigenvalues. Thus for example in three dimensions, we have

λ = 2 1 + √ d x 2 1 + d y 2 1 + √ d z . ( 30
)
The size of this eigenvalue is the established measure for the content of condensate in the wave function. The remaining eigenvalues can be found in terms of the largest one:

λ n = λ (1 -λ) n , ( 31
)
where n is a non-negative integer [14]. We thus see that if λ is close to one then all other eigenvalues will fall off very fast, whereas for small λ one has a distribution of many non-zero eigenvalues and a highly fragmented state.

Symmetries imposed by boson and fermion statistics

The solutions discussed above are all obtained without any requirements of symmetry due to groups of identical particles in the system. The hamiltonian necessarily must commute with correponding permutation operators, and the solutions should then have the proper symmetry. It is then only a question of selecting those states among the complete set of all solutions. Unfortunately, this tempting conclusion is wrong in general because it is based on the assumption of non-degenerate states. Degeneracies allow mixing of different symmetries, which implies that the solutions are linear combinations of the required symmetries. Thus, a simple method to restore the required symmetry is to construct linear combinations of the available solutions with the same energy.

We have assumed that the hamiltonian is independent of intrinsic properties like spin of the particles. This means that the total wavefunction factorizes into intrinsic (including particle spins) and spatial parts. Here we only consider the symmetries of the spatial part which subsequently has to be combined with the remaining parts of the wavefunction. The task is not easily formulated in general since the system may consist of different groups of particles each with their own symmetry requirement. This could be a combination of identical bosons and fermions, or identical particles placed in different geometries by external fields and consequently effectively behaving as nonidentical particles.

To illustrate the method we assume that all particles are identical, or alternatively we omit reference to the non-identical particles left in their orbits while only considering those explicitely mentioned in the following formulation. We write the symmetrized or antisymmetrized spatial wavefunction Φ

Φ( r 1 , r 2 , ..., r N ) = N p p π p Ψ( r p(1) , r p(2) , ..., r p(N ) ) , (32)
where r i = (x i , y i , z i ) is the coordinate for particle i, p is a permutation of the set of numbers {1, 2, ..., N}, N p is a normalization constant and π p is +1 for spatial (boson) symmetry, and π p = ±1 for even and odd permutations p, respectively, when we construct wavefunctions with spatial (fermion) antisymmetry. The energy of Φ is given as the energy, E n x,k ,n y,k ,n z,k from Eq.( 21), which is the same for each term of Φ. This conclusion is due to the fact that the hamiltonian is invariant under these coordinate changes, and the energy of the state in Eq.( 32) is the expectation value of the hamiltonian. The wavefunction in Eq.( 32) may be identically vanishing, e.g. when a symmetric, Ψ, is antisymmetrized or vice versa, when an antisymmetric, Ψ, is symmetrized. In these cases the state, Φ, does not describe any physical system. The procedure can be specified in more details with matrix manipulations. If the coordinate transformation corresponding to one of the permutations in Eq.( 32) for one of the Cartesian coordinates, x, is described by P x , i.e. x p = P x x (with x = (x 1 , x 2 , ..., x N ) T ), we can express the equivalent transformation in the final coordinates by v p = MP x M -1 v, with M defined in Eq. (19). Thus by replacing v by v p in the arguments of the terms in the wavefunctions in Eq.( 32) the same variables v are used as in the unpermuted term. The same coordinates are then used to express all terms in the full (anti)symmetrized wavefunction

The permutations leave the exponent in the wavefunction unchanged, i.e. v T p B p v p = v T B v. The reason is that the N -1 oscillator frequencies, masses and lengths are identical. The last length is related to the center of mass motion, which either has to be omitted without confining field or remains invariant under the permutations as decoupled from the intrinsic motion described by the N -1 frequencies. The Hermite polynomials change arguments but remain polynomials of the same order in the new coordinates v p . In total, each term in Eq.( 32) is a different polynomial of the same order in the coordinates corresponding to the permutation. The expectation value of an arbitrary operator is then very difficult to find due to the many terms. Operator expectations are then analytic if they can be found in an oscillator basis, and in particular this is the case for any polynomial structure of the coordinate or momentum operators. However, in general the expectation values may consist of many terms.

Small excitations leave only few different oscillator quanta in the wavefunction. The non-vanishing terms in Eq.( 32) are then limited to rather few. The extreme case is the state with all n x,k = 0 which is completely symmetric under all permutations, that is only one term describing the ground state for identical bosons. This case will be our main concern in the remaining part of the paper. When all quanta are equal to zero except one with an arbitrary n x,k = 0, the number of different permutations are N , which is a manageable number.

Another practical limit is when the number of identical particles effectively is small, e.g. when identical particles are confined by external fields to different states or spatial locations. Then only very few permutations are present, and the brute force method is easily applicable. In any case the brute force methods are only necessary when information beyond the energy is required. More suitable procedures can no doubt be designed for specific systems and corresponding operators, as for instance discussed recently in the context of the hyperspherical harmonic approach [15]. However, all derivations can still be made analytically if the operator expectations are calculable for oscillator states.

Bosons interacting in a trap

Two particles interacting in a trap is obviously the simplest non-trivial system of N interacting particles. Still very interesting features for two particles were discovered for the extreme limit of a contact interaction and a confining trap potential [9]. It would not be surprising if an oscillator approximation has difficulties in a quantitative description of such systems. On the other hand this is a challenging problem, and the structure of N particles in a trap with such pairwise interactions are an active field of research. We shall therefore attempt to extract systematic overall features within our analytical formalism. The philosophy is to reproduce the two-body properties as close as possible, and then systematically calculate properties of the many-body system.

Several other approaches have been used for few-body systems of atoms interacting via the contact interaction. Among these are effective field theory approaches [16,17,18,19], shell-model [20,21], and Monte Carlo calculations [22,23,24] in three dimensions for fermions, various hyperspherical and variational approaches in three [25,26,27,28,29] and in two dimensions for bosons [30,31,32,33], and exact diagonalization in two dimensions for fermions [34,35] and for bosons [35]. The zerorange interaction does not allow the diagonalization of the many-body Hamiltonian in spatial dimensions greater than one, so all these methods must address that in some way. Most often, this manifests as a cut-off to a particular subspace, which effectively renormalizes the interaction strength [16,36,20]. In [33], a Gaussian form of the interaction is used to approximate the zero-range potential, and only repulsive interactions are considered. With the exception of [34], all discuss energy results as a function of model space size, and do not yet discuss other observables.

Adjusting the oscillator parameters

First we must decide how to choose the parameters in the oscillator model. For two identical bosons we have initially five parameters, i.e. mass, interaction frequency, energy shift, confining frequency, and center position of the confining field. One of these can always be chosen as a scale parameter or a unit without consequences for any of the properties. In the present case, the external field is provided by its frequency, ω 0 , independent of the two-body interaction, which then leave the interaction frequency, ω ik = ω for all i, k, and the energy shift, V ik = V for all i, k, to be determined. These parameters all identical since we consider indistinguishable particles. We can also immediately set the center position of the confining field to zero, as there is no such center shifting or multi-center effect in the external field in the original two-body problem.

Our strategy will be to reproduce the ground state properties of the model of Busch et al. [9] as much as possible using a harmonic oscillator. We work exclusively in the domain where the lowest state is the molecular branch that represents a bound state when the trap is removed. The population of this particular state was considered previously in [37,38]. In three dimensions, this means that we consider only the positive scattering length side of the resonance, whereas in two dimensions this branch is always present. This procedure provides our effective mapping from the two-body interaction to the solvable N -body problem. We note that our choice of branch means that for large scattering lengths the two-body energy goes to 1/2 ω 0 (three dimensions) and ω 0 (two dimensions), where ω 0 is the external field frequency. For small scattering lengths it goes as the inverse of the scattering length squared in both cases, since it represent the universal bound state that is also present when ω 0 → 0.

The Hamiltonian, solved in [9] for two particles, is

H = - 2 2m ∇ 2 r + 1 2 mω 2 0 r 2 + 4π 2 a m δ (3) (r) ∂ ∂r r , ( 33
)
where m is the mass of the particles, ω 0 is the external trap frequency as mentioned above, r is the relative coordinate r = 1/2(r 1r 2 ) and a is the scattering length of the two-body potential assumed to be a regularized δ-function.

The solutions are given as eigenvalue equations and corresponding wave functions, ψ. The form of ψ for both two and three dimensions is found to be

ψ(r) ∝ 1 2 π -D/2 e -r 2 /(2 2 ) Γ(-ν)U (-ν, D/2, r 2 / 2 ) , (34)
where the dimension is D = 2, 3, the length is given by 2 = 2 /(mω 0 ), and the relative energy E rel /( ω 0 ) = 2ν + D/2 is given in terms of the non-integer quantum number, ν.

The eigenvalue equations are respectively

2 Γ(-ν) Γ(-ν -1/2) = a , D = 3 ( 3 5 ) γ + 1 2 ψ(-ν) = ln a , D = 2 , ( 36
)
where γ is the Euler-Mascheroni constant. For a given scattering length and trap frequency ν is obtained and both energies and wave functions are determined. There are many solutions to the above equations, but we repeat that we work exclusively with the lowest molecular bound state, corresponding to the lowest solution for ν. Pertinent features from these solutions are now used to choose the oscillator parameters. First we directly choose the same external frequency ω 0 for both particles. Second, we compute the mean square radii for D = 2, 3 from (34) and equate to the corresponding oscillators, i.e.

ψ|r 2 |ψ ψ|ψ = D 2µ ω 2 + ω 2 0 . ( 37
)
This determines the oscillator frequency ω. Finally, we adjust the energy shift for the oscillator model to reproduce the correct two-body energies, i.e.

(2ν + D/2) ω 0 = D 2 ω 2 + ω 2 0 + V , ( 38
)
where ν is obtained by solving the relevant eigenvalue equation from (35) and (36). The energy shift is V shift (N = 2) = V as seen in ( 5).

The interaction frequency combined with the trap frequency determines all structure. The size of the system is crucial and we have chosen to reproduce the radius. It is not as meaningful to adjust to energies since the oscillator cannot be expected to describe very weakly bound and spatially extended structures. Still we expect to get an indication of the energy variation with N from the shift in the energy zero point. All oscillator parameters are now determined for identical bosons, and we can proceed to investigate consequences for the many-body system.

We note the method to determine the oscillator parameters for the two-body potential can be considered much more general and other states may be used, like for instance a higher excited state in the zero-range model. In that case certain technical considerations arise, for example with respect to the nodal surfaces of the many-body wavefunction in relation to the two-body interactions from which it is built. This is particularly important in the case of fermions due to the requirement of antisymmetry. A discussion of these questions will be presented elsewhere.

The N -body system

The total energy shift for N particles is simply the number of pairs times the two-body shift, i.e.

V shift = 1 2 N (N -1)V . (39)
The external frequency is ω 0 for all particles. Solving the oscillator model leads to a set of frequencies where one of them equals the external trap frequency and corresponds to the decoupled center-of-mass motion. The remaining N -1 other frequencies are degenerate and for each direction given by ωx = N/2 ω 2 + 2ω 2 0 /N .

The total ground state energy of the relative motion is then

E gs = 1 2 N (N -1)V + D 2 (N -1)ω x . (41)
For repulsive interactions the sign of ω 2 in ( 40) and ( 41) is changed. If ωx becomes imaginary the system becomes unstable, i.e. if Nω 2 < -2ω 2 0 , then the pairwise repulsion is too strong for the external field to confine the system.

The energy expression in (41) for a gas of many identical bosons has N -dependent terms of different origin. The term proportional to N 2 is solely from the interaction, whereas the N 3/2 term originates from kinetic energy, two-body interaction, and external one-body potential, as given by the solution. The relative influence of the external potential decreases with N . The remaining part still varies as N 3/2 which is a compromise between the pairwise interaction increasing as N 2 and the linear kinetic energy depending on N .

We can compare with energy relations directly derived for a gas of bosons that interact via a δ-function two-body potential [39,40], which in the limit of weak interaction in three dimensions is

E 3N 2 ω 0 + N 2 U 0 2(2π) 3/2 b 3 , U 0 = 4π 2 a m , ( 42
)
where b is the external trap length. This expression has the same N 2 scaling from the interaction as V shift in the oscillator model. The proportionality factor in V shift is obtained through the frequency and strongly depends on which state is used for the adjustment. A more detailed comparison is then less direct. The linear term in ( 42) is obviously arising from kinetic energy and external field which for the oscillator leads to the term proportional to N 3/2 . For a strongly interacting system, where Na/b 1, a variational calculation where part of the kinetic energy is neglected gives a different dependence [39,40]

E = 5 4 2 π 1/5 Na b 2/5 N ω 0 . ( 43
)
The overall energy scales as N 7/5 which presumably should be compared to our result when ω ω 0 where we get a linear energy scaling from kinetic energy and part of the two-body interaction. Another part of the interaction is still scaling as N 2 . A similar result is not surprisingly found with Thomas-Fermi approximation, where the kinetic energy term is ignored, but where higher order terms can be included to improve the description [41,42,43].

The mean square distance from the center-of-mass is also calculated to be

(r -R) 2 = D(N -1) 2N 3/2 m 1 ω 2 /2 + ω 2 0 /N , ( 44
)
where R is the location of the center-of-mass. This measure of the spatial extension is a reflection of the inverse behaviour of energies and mean square radii.

Energies and radii

The dependence on particle number N is very explicit, but still two or rather three terms compete with their different N -scaling. We first show the relative energy per particle (i.e., the energy of the relative motion of the particles without the center-ofmass contribution corresponding to one particle moving in the external field) in figure 1 as function of N for three dimensions. It happens that a two-body system is "unbound" in the oscillator model corresponding to positive relative energy. This occurs when the positive contribution in (38) is larger than the negative V shift . However, V shift must dominate as N increases, and in fact we find that once another particle is added, the relative energy is less than zero for any of the studied scattering lengths. The N -dependence is very smooth and steadily increases the binding which varies strongly when the scattering length approaches zero in units of the trap length. This is also the region where the magnitude of the two-body bound state energy increases rapidly. The frequency dependent term in the energy proportional to ωx is simply the zero point motion of an oscillator. It is therefore also the smallest unit of excitation of the system corresponding to one particle lifted in one dimension from the ground to first excited state. This is the energy of the normal mode of excitation. The dependence is shown in figure 2 as function of scattering length from the Busch et al. model in (33). The frequency is small and constant for scattering lengths larger than the trap length, while it begins to grow very quickly as soon as the trap length exceeds the scattering length.

The size of the system is, along with the energy, the most fundamental property. The intuitive implication that smaller radii follow larger binding is also observed in figure 3, regardless of the size of the scattering length, though at the weaker interactions the change is rather flat for the first few added particles. The N -dependence is again rather simple as seen in (44) where the mean square radius decrease with 1/

√

N for large N . Otherwise the radii are varying with frequency as usual. The interesting part is rather that the sizes increase substantially with scattering length for given particle number.

We now repeat the procedure in two dimensions. The major difference is obviously from the adjusted oscillator input parameters. The energy and radius expressions are already given in (41) and (44). The N -dependence of the energy is shown in figure 4 where the same overall structure as in figure 1 appear. As in 3D, the energy decreases monotonically with the addition of more atoms for all scattering lengths.

The two-body energy shift is less negative because the size requirement to determine the frequency better matches the one-body frequency. Then several of the N -body scattering length as a function of particle number. From bottom to top, the ratios are a/ = 100, 10, 5, 2, 4/3, 1, 1/2, 1/5, and 1/10. Note that a few of the bottom lines at large scattering length over trap length ratios appear to stop at small particle numbers. This is because there is a sign change in the energy, which is discussed around (45).

energies are positive, although "binding" (negative energy) finally occurs by adding a sufficient number of particles. This critical number, N c , where the system becomes self bound is given by

N c = ( ω) 2 V 2   1 + 1 + 4 ω 2 0 ω 2 V ω 2   , (45)
which for the two dimensional case depends strongly on the initial scattering length as seen in figure 5. For the largest scattering length in figure 5 N c is about 12, whereas there is binding for all particles for scattering lengths smaller than about three. As the scattering length becomes arbitrarily large, the critical number also increases, meaning that the addition of any number of particles is not sufficient to bind the system at infinite scattering length.

In figure 6 we show how the degenerate frequency for ten particle behaves as a function of the scattering length. The behaviour is very similar to then 3D result, with both curves turning up strong after the scattering length becomes less than the external confinement length.

Figure 7 shows the size results in two dimensions. These also follow a slightly different trend than in three dimensions. In two dimensions, for several scattering lengths, the size actually increases for the addition of a particle (going from three to four particles), before falling for all additional particles. The scattering lengths where this behaviour is seen correspond roughly to those that have positive three-body energies. For smaller scattering lengths, the relative sizes decreases monotonically with an increase in particle number.

One-body density

The one-body density matrix has information about the mean-field content of the wave function [11]. This is quantified by the eigenvalue in (30), which directly measures how much this state has the structure of a coherent (condensed) state. In figure 8, we see how this eigenvalue, λ, evolves for 3D with interaction strength and particle number. For the most part it increases with particle number, though at weaker interactions there is a small minimum around five or six particles and it increases thereafter, while it decreases uniformly with interaction strength. The overall increase with N is in agreement with the theorem that the mean-field wave function is approached as N tends to infinity [11]. The condensate fraction is large for large scattering lengths where the external field is decisive, and consequently favors the corresponding mean-field structure. This follows from the fact noted above that the two-body wave functions become essentially non-interacting oscillator states given by the confinement (ω becomes small). On the other hand, for small scattering lengths the structure is far from that of a condensate. The particles are much more tightly bound with strong correlations. Again this is consistent with the fact that we have very strongly bound two-body states in this limit that are very different from the non-interacting harmonic states.

The question of whether a true condensate actually exists in a realistic quasi-2D setup with harmonic trapping potentials below a certain critical temperature [44,[START_REF] Bloch | [END_REF] will not be addressed here. We simply take the appearance of a large eigenvalue in the one-body density matrix as our working definition of a condensate as in the 3D case above. The condensate fraction are shown in figure 9 for the 2D system. The same tendency of increase with N is present here. Again the condensate fraction is large (small) for large (small) scattering lengths compared to the size of the external field. This reflects the amount of correlation in the wave function precisely as for 3D.

Quantitatively we find that λ is even flatter as a function of particle number than in 3D, and it is also consistently higher for the same scattering length. We note that the 2D results consistently produce larger condensate fractions than in 3D. As we discuss below this is connected with the fact that the 2D case resembles a non-interacting system when a → ∞ better than the 3D case does. The remaining deviation of λ from unity is due to the separation of center-of-mass.

Normal modes and symmetries

The characteristic properties of a system are reflected in the structures and energies of the normal modes. Here we discuss the energies of the degenerate frequencies which is the amount of energy required to excited the modes. Unfortunately, the degeneracy is itself an obstacle for understanding the uncoupled structures of single excitations. The reason is that any wave function expressed as a linear combination of the same energies is equally well qualified. The set of normal modes only has to be orthogonal, but that can be achieved in infinitely many ways. One general exception is if other conserved quantum numbers have to be restored by specific linear combinations of the degenerate states. Angular momentum is a prominent example in the absence of external fields.

We now discuss what determines the degeneracy, or equivalently what would break it. First, if all masses, all two-body interaction frequencies, and all one-body frequencies are equal, then N -1 degenerate frequencies are produced along with the one-body frequency, which is returned unchanged. If masses are changed, then the degeneracies are broken, and if N -2 masses are changed then all degeneracies will be broken. This is the same for the one-body frequencies (those in (3)), if N -2 of them are different from each other, then all the output frequencies will be different.

For the interaction frequencies in (2) the situation is slightly more complex. If one row and column in the A-matrix of (6) have the same interaction frequency, i.e., if this particle interacts the same with all other particles, then there will be at least one pair of degenerate frequencies. In general, of the N (N -1)/2 interaction frequencies (off-diagonal elements of the A-matrix), if (N -1)(N -4)/2 (N > 4) of them are different, then that is enough to guarantee all symmetries are destroyed. However, all degeneracies can be broken with a wiser distribution of the different frequencies. If, for example, the frequencies immediately above the main diagonal of A are all different and all the remaining ones being the same (a total of N different off-diagonal frequencies), then that is enough to destroy all degeneracies of the resulting normal modes.

Thus, degeneracy can be broken or reached through many different paths. The structure of the resulting degenerate normal modes depends on the chosen path. We still find it interesting to investigate specially selected normal modes. If all particles are distinguishable it should in principle be possible to observe corresponding vibrational structures where each particle is detected. To probe the underlying structure, revealed by the normal modes, we therefore approach the degenerate limit from an entirely nondegenerate system. The two simplest ways are to differentiate the particles by minute differences in mass, or alternatively in one-body frequency.

We first notice that the normal modes are the results of transforming the set of oscillators to diagonal form. The energies and eigenvalues are obtained from the matrix depending on masses, external field frequency, and two-body interaction frequency.

Thus by choosing all these parameters to be identical the normal modes in two and three dimensions are the same. This is achieved by the same masses, same external field, and corresponding choices of scattering lengths in 2D and 3D such that the twobody interaction frequencies from (40) are identical. The correspondence is seen in figure 10 where both scattering lengths are small at the same time and also increase simultaneously. The 2D scattering length approach an upper limit of about 1.14 2D for increasing a 3D . The size in 2D is much smaller than in 3D. In figure 10 we also show the relation between 2D and 3D oscillator shifts obtained from (38) for corresponding scattering lengths. It is clear from the figure that it is possible to map 3D results onto 2D results, but the reverse is not always true. This is due to the different behaviour of the systems for large scattering length, where the energy and square radius of the two-body system are controlled by the properties of the external trap. In two dimensions we find that r 2 = 2 when a → ∞ and (37) then tells us that we have to choose ω = 0, i.e. we have a non-interacting system. For three dimensions, the situation is somewhat different as the virial theorem applied for a → ∞, tells us that r 2 = 2 /2. From (37) we deduce that ω = 2 √ 2ω 0 , i.e. our harmonic equivalent is an interacting system still. This also implies that the interaction frequency in 3D, ω in (40), can never be smaller than 2 √ 2ω 0 . In 2D, there is no lower limit for the interaction frequency and it eventually vanishes at a large enough scattering length where the external field determines all properties.

From figure 10 we can in principle from a series of 2D calculations for different scattering lengths extract 3D results. The procedure is to start with the desired 3D scattering length, find the corresponding 2D scattering length, read off the related 2D results for normal modes, radii, and oscillator energy shift. The normal modes and radii are now the 3D results, and the 3D energy is obtained by finding the related 3D oscillator energy shift from figure 10 combined with the expression in (41). The procedure can be reversed provided the desired 2D scattering length is within the range accessible to conversion from 3D (i.e., a 2D / 2D ≤ 1.14).

The hamiltonian we obtain in our oscillator approximation can be decoupled by a change of coordinates as discussed above. The normal modes can therefore be viewed in one dimension at a time as the amplitudes with which each of the individual particles are moved when the corresponding mode is excited on top of the ground state. More explicitly, consider exciting the i'th mode with probability p, so that the wavefunction of the system becomes |Ψ(t) = 1 -p 2 |0 + √ pe -iω i t |i , where ωi is the frequency of mode i. The displacement in time is then x i (t) = Ψ(t)|x|Ψ(t) = A 0i cos(ω i t), with amplitude given by A 0i = 2 p(1 -p 2 ) 0|x|i . We illustrate the modes pictorially below. The fact that the different spatial directions are decoupled also implies that a spherical external field produces degenerate modes. This degeneracy can be lifted by deforming the field but such a symmetry breaking would only reflect the properties of the field.

We therefore consider the one-dimensional eigenmodes which apply for both 2D and 3D when using a correspondence of scattering lengths down to 1D as was done between 2D and 3D in figure 10.

In figure 11 and figure 12 we show a sequence of modes for six and eight particles respectively. The dependence on scattering length is rather weak for both small and large a/l. The general picture that emerges from breaking the symmetry by mass differences is first that the energy of the center-of-mass oscillation is maintained. Second, the energies of the smallest and the largest mode correspond to oscillation where either the lightest or the heaviest particles move in one direction while all the others move less and in the opposite direction. The remaining modes with intermediate energies correspond to one particle joining the lone particle, moving in the order of increasing mass for each mode, though in most cases it appears that only one or two particles are displaced significantly, while the remainder stay in a slightly displaced clump. It should be noted that the normal modes are highly sensitive to the symmetries of the interaction, and how we broke the symmetry in order to show non-degenerate normal modes. In this present case, we changed the masses slightly to make the particles distinguishable, but still with identical interactions between all particles. We expect that if the interactions are changed slightly in a certain prescribed manner, then the normal modes will reflect the symmetries of that change.

Summary and Conclusion

The properties of the quantum mechanical N -body system are determined from the basic one-and two-body interactions. However, in general this problem is very hard Amplitudes of the normal modes in the x direction of a six particle system in 3D calculated with a/ = 1, which is equivalent to a ratio of 0.61 in the 2D. The largest amplitude is rescaled to unit magnitude. Points in the diagram are numbered from one to six with the mass increasing from particles one to six. Overlapping points are not labeled, and appear in groups usually close to the origin, which is indicate by the black solid horizontal line. The points oscillate through the equilibrium position (the black solid line) with a frequency of the given normal mode(e.g., x i (t) = A 0i cos(ω i t), where A 0i are the amplitudes shown in the figure). The frequencies from left to right are 1.003, 10.433, 10.447, 10.456, 10.470, and 10.487 in units of ω 0 and are not plotted to scale. The center of mass mode is the first mode shown. Amplitudes of the modes in the x direction of an eight particle system in 3D with a/ = 1, which has an equivalent ratio in 2D of 0.61. As in figure 11, the largest amplitude is normalized to one unit of length. Points in the diagram are numbered from one to eight with the mass increasing from particles one to eight (overlapping points are not labeled). These points oscillate through the equilibrium position, indicated by the solid black horizontal line, with a frequency of the given normal mode. The frequencies of the normal modes are, from left to right: 1.004, 12.030, 12.040, 12.050, 12.059, 12.069, 12.085, 12.097 in units of ω 0 and are not plotted to scale. The center of mass mode is the first mode shown.

or impossible to solve. Here we use an approximate approach which replaces the interactions with quadratic forms in the coordinates, either by direct fits of the potentials or by adjusting parameters to reproduce crucial properties. This approximation allows analytical investigations of the N -body system with the properties expressed in terms of the two-body characteristics. Having such a direct approach to the general many-body problem can provide both important analytical insight and be a valuable benchmark for more intricate methods.

The most general harmonic oscillator potentials to describe one-and two-body interactions allow analytical solutions of the N -body Schrödinger equation. However, the center-of-mass degree of freedom requires special attention as we have described above. We have employed Cartesian coordinates and one-, two-, and three-dimensional systems are therefore equally simple to handle. We have discussed the properties of the resulting set of decoupled oscillators and their relation to the initial interactions. More explicitly we calculated energies, radii, and the one-body density matrix and presented its spectrum.

As an application of our formalism, we consider first the simplest case of N bosons in a trap interacting through contact potentials. This was done by using a mapping from the energies and radii for the two-body system to the oscillator parameters for the analytical calculations. We note that while our choice of energy and radius as the essential parameters to reproduce in the two-body system was perhaps the most physically reasonable one, other choices are also possible as long as one has two quantities that fix the oscillator frequency and the shift. We calculated energies and radii as function of boson number and the scattering length as a measure of the interaction strength. Typically the binding energies increase and the radii decrease with N . We also discussed characteristic limits for scattering lengths much smaller and much larger than the length scale of the external trap.

An interesting calculation was done for the one-body density matrix for which we calculated the dependence of the lowest eigenvalue with the number of bosons and the scattering length, with a careful treatment removing the center-of-mass degree of freedom. The limits of large scattering length gave large condensate fractions, whereas for small scattering length we found a fragmented state. The former is due to the near non-interacting system, whereas the latter is caused by the strongly bound molecular state which introduces substantial amounts of correlations in the system. These conclusions apply equally for both two and three dimensions. We also computed the critical number of particles that can form a self bound system of negative energy. This number increases with increasing scattering length. In three dimensions this happens already for a few particles while in two dimensions more than 50 bosons are needed for large scattering lengths.

As a novelty, we considered the normal modes which are characteristic properties of any system. However, in the case of the N -boson system, there is a large degree of degeneracy that can obscure the detailed behaviour. The ambiguity due to degenerate eigenmodes is circumvented by breaking the symmetry and then approaching the limit of full degeneracy for the identical boson system. A different way of breaking the symmetry is to deform the external field but the resulting eigenmodes then reflect precisely the chosen deviation from spherical symmetry. We then computed the one-dimensional oscillatory eigenmodes which can be related to equivalent values of the scattering lengths in two and three dimensions. What we found in the normal modes was an interesting tendency for the particles to cluster in smaller groups and then perform motion with respect to other such groups. This implies that excitations can induce strong correlations in a many-boson system even when all interactions are equal.

The method that we have presented in this report is completely general and can be applied to systems in external fields or to self-bound structures. The treatment of deformation of the trap or the two-body interaction is straightforward and we expect that rotation of the external trapping potential poses no problem as well. While we have only treated bosons in this work, the extension to fermions is achieved through proper (anti-) symmetrization of the wave function and also Bose-Fermi mixtures are accessible. With the possibility of having displaced centers we can also apply the method for split traps, and also even more exotic geometries with mixed dimensions which are under current study within Fermi-Fermi mixtures of ultracold atoms [45,46,47,48]. Our method can also be applied to cold polar molecules, in particular to the case where two-dimensional confinement is induced to make layered systems [49] where non-trivial two-and three-body bound states appear in the bilayer [50,51,52,53,54] that open for study of various exotic many-body states in both bilayer and multi-layer systems [55,56,57,58,59,60]. The form of the dipolar potential has harmonic oscillator shape in the inner part when particles are confined in multi-layers and we therefore expect the the method presented here to readily provide analytical results valid for intermediate and strong dipolar strength. The fact that normal modes are easily accessible with the harmonic method presented makes this even more interesting. This can help understand the modes of excitation in chains and complexes for use in thermodynamic calculations of system properties.

In conclusion, the analytic solutions of coupled harmonic oscillators can be used to study the overall properties of N -body systems and structures can be calculated in general from two-body properties for different particle number and geometries. This can serve as a valuable complement to the understand of results obtained with more intricate methods or help solve systems that are intractable in other approaches.

Figure 1 .

 1 Figure 1. Relative energies per particle at several different values of the 3D scattering length as a function of particle number. The different values of the ratio are, from bottom to top: a/ = 100, 10, 5, 2, 4/3, 1, 2/3, 1/2, 1/5, and 1/10.

Figure 2 .Figure 3 .

 23 Figure 2. Degenerate frequency, ω of the N -particle systems in 3D plotted logarithmically as a function of the ratio of the scattering length to the external confinement length. The different values of N plotted are, from bottom to top: 3, 4, 5, 6, 10, 20, and 30.

Figure 4 .

 4 Figure 4. Relative energies per particle at several different values of the 2Dscattering length as a function of particle number. From bottom to top, the ratios are a/ = 100, 10, 5, 2, 4/3, 1, 1/2, 1/5, and 1/10. Note that a few of the bottom lines at large scattering length over trap length ratios appear to stop at small particle numbers. This is because there is a sign change in the energy, which is discussed around(45).

Figure 5 .Figure 6 .

 56 Figure 5. The critical number N c for 2D as function of the ratio of the scattering length to the external length. The vertical line is a/ = 3.3, at which point (to the left of this line) all systems of more than two particles become self-bound.

Figure 7 .

 7 Figure 7. Relative sizes at several different values of the 2D scattering length as a function of particle number. From top to bottom, the ratios plotted are a/ = 100, 10, 5, 2, 4/3, 1, 1/2, 1/5, 1/10.

Figure 8 .

 8 Figure 8. The value of λ at several different values of the 3D scattering length as a function of particle number. From top to bottom, the ratios plotted are: a/ = 100, 10, 5, 2, 4/3, 1, 2/3, 1/2, 1/5, 1/10.

Figure 9 .

 9 Figure 9. Condensate fraction at several different values of the 2D scattering length as a function of particle number. From top to bottom, the ratios plotted are: a/ = 100, 10, 5, 2, 4/3, 1, 1/2, 1/5, 1/10.

Figure 10 .

 10 Figure10. The relation between the 2D scattering length to external length ratio to the same quantity in 3D for both the interaction frequency (red curve) and energy shift (green curve).

Figure 11 .

 11 Figure11. Amplitudes of the normal modes in the x direction of a six particle system in 3D calculated with a/ = 1, which is equivalent to a ratio of 0.61 in the 2D. The largest amplitude is rescaled to unit magnitude. Points in the diagram are numbered from one to six with the mass increasing from particles one to six. Overlapping points are not labeled, and appear in groups usually close to the origin, which is indicate by the black solid horizontal line. The points oscillate through the equilibrium position (the black solid line) with a frequency of the given normal mode(e.g., x i (t) = A 0i cos(ω i t), where A 0i are the amplitudes shown in the figure). The frequencies from left to right are 1.003, 10.433, 10.447, 10.456, 10.470, and 10.487 in units of ω 0 and are not plotted to scale. The center of mass mode is the first mode shown.

Figure 12 .

 12 Figure12. Amplitudes of the modes in the x direction of an eight particle system in 3D with a/ = 1, which has an equivalent ratio in 2D of 0.61. As in figure11, the largest amplitude is normalized to one unit of length. Points in the diagram are numbered from one to eight with the mass increasing from particles one to eight (overlapping points are not labeled). These points oscillate through the equilibrium position, indicated by the solid black horizontal line, with a frequency of the given normal mode. The frequencies of the normal modes are, from left to right: 1.004, 12.030, 12.040, 12.050, 12.059, 12.069, 12.085, 12.097 in units of ω 0 and are not plotted to scale. The center of mass mode is the first mode shown.


R Wiegner
email: ralph.wiegner@physik.uni-erlangen.de

J Von Zanthier

G S Agarwal

Quantum interference and non-locality of independent photons from disparate sources

Keywords: numbers: 42, 50, Dv, 03, 65, Ud

We quantitatively investigate the non-classicality and non-locality of a whole new class of mixed disparate quantum and semiquantum photon sources at the quantum-classical boundary. The latter include photon added thermal and photon added coherent sources, experimentally investigated recently by Zavatta et al. [Phys. Rev. Lett. 103, 140406 (2009)]. The key quantity in our investigations is the visibility of the corresponding photon-photon correlation function. We present explicit results on the violations of the Cauchy-Schwarz inequality -which is a measure of nonclassicalityas well as of Bell-type inequalities.

I. INTRODUCTION

The question of interference between independent photons has attracted considerable attention since the celebrated work of Hong, Ou and Mandel who showed how two photons of the pair emitted in a non-linear crystal by the process of spontaneous parametric down conversion (SPDC) can interfere on a beam splitter [1]. Many subsequent investigators verified the Hong-Ou-Mandel interference using a variety of SPDC sources [2][3][4]. Also the interference of photons from truly independent quantum sources, e.g., two distant atoms, has been analyzed in great detail [5,6] and recently observed [7][8][9][10]. Note that the interference is displayed in terms of the photonphoton correlation function [11] rather than in measurements of the mean intensity. So far, most works have concentrated on the interference produced by independent but identical sources. However, the questions then arise: what is the extent to which photons from independent but disparate sources can interfere [12]? Further, what is the nature of the radiation field generated and to what extent are the properties of the field strictly quantum? Can one use such a quantum character to shed light on the non-local character of the field? We provide quantitative answers to these questions. Specifically we investigate the quantum interference of two disparate sources, one being a single photon source like an excited atom and the other a source which can give rise to properties ranging from completely classical to completely quantum. In the latter category we specifically examine sources which are either in a state called photon added coherent state [13] or photon added thermal state [14]. The quantum properties of these sources arise from the fact that one has added photons to classical states like thermal and coherent states. Such sources have been recently experimentally realized and their properties studied [15][16][17][18]. We calculate the interference pattern of the cor-responding photon-photon correlation functions, derive their visibilities and discuss the conditions under which the photon-photon correlation functions displays strictly quantum behavior. The latter is formulated quantitatively in terms of the Cauchy-Schwarz (CS) inequality [19] valid for classical fields. After revealing the quantum nature of the light fields generated by disparate sources we discuss the conditions under which Bell inequalities can be violated [20][START_REF] Bell | Speakable and Unspeakable in Quantum Mechanics[END_REF][START_REF] Clauser | [END_REF][23].

The organization of the paper is as follows: In Sect. II we derive the photon-photon correlation functions and the corresponding visibilities for combinations of mixed disparate quantum-classical and quantum-semiquantum light sources. As quantum source we consider throughout the paper a single photon source, e.g., an initially excited two-level atom spontaneously emitting a single photon. As second source we assume either a purely classical source exhibiting Poissonian or thermal statistics [24][25][26] or a semiquantum source emitting a field which is either in a photon added coherent state (PACS) or a photon added thermal state (PATS). In Sect. III we formulate a version of the CS inequality for position dependent photon-photon correlation functions in order to investigate quantitatively when the radiation fields display non-classicality, by explicitly showing the regimes where the CS inequality is violated. In Sect. IV we then discuss, under which conditions Bell inequalities can be violated by our mixed quantum-classical and quantumsemiquantum sources. In Sect. V we conclude with a summary of our results. two detectors, i.e., we assume that each of the two detectors registers one photon so that two photon absorption processes at one detector can be excluded. The detectors are located at positions r 1 and r 2 , where |r i -R n | d (i = 1, 2, n = A, B), so that the far field condition is fulfilled. To simplify further calculations we only consider coincident detections, i.e., simultaneous events at both detectors.

FIG. 1: Setup of the investigated arrangement. We consider one photon source localized at position RA and another photon source at position RB. Using two detectors at positions r1 and r2 we assume that each detector registers a photon in the far field of the sources.

In Glauber notation [11] the photon-photon correlation function can be written as

G (2) (r 1 , r 2 ) = Ê(-) A+B (r 1 ) Ê(-) A+B (r 2 ) Ê(+) A+B (r 2 ) Ê(+) A+B (r 1 ) , (1)
where ( Ê(-) A+B ) † = Ê(+) A+B and . . . denotes the expectation value corresponding to a temporal average for stationary fields. Hereby, the positive frequency part of the electric field operator Ê(+) A+B (r i ) takes the form [27] Ê(+)

A+B (r i ) ≡ Â(+) i + B(+) i ∝ e -i(kri RA+ψA) s + e -i(kri RB+ψB ) a, (2)
where k = 2π λ = ω c denotes the wave number of the two sources, assumed to radiate at the same frequency ω, and ri := ri |ri| is a unit vector in the direction of the ith detector. ψ A and ψ B are the initial phases of the independent sources A and B, which in general may be time dependent.

Except for the next paragraph, we consider in this paper for source A a single photon source, e.g., a two-level atom with upper level |e and ground state |g . The atom is assumed to be initially fully excited to the state |e so that after the spontaneous emission of a single photon the atom is transferred to the state |g . In this case the operator s in Eq. ( 2) denotes the lowering operator |g e|. As source B we consider in the following either a more traditional source like a coherent source or a thermal source [24][25][26] or a source which displays both, classical and quantum properties. For the latter we specifically investigate the configurations where source B is either in a PACS [13] or in a PATS [14]. The operator a in Eq.

(2) thus denotes the photon annihilation operator for the source B.

For a better comparison with the following calculations we start to review the pure classical case where the two independent light sources A and B emit coherent light [5]. The operator s in Eq. ( 2) denotes in this case the photon annihilation operator for the source A. Since the two light fields are in a coherent state, |α A and |α B with mean photon numbers nA and nB , respectively, we calculate from Eqs. ( 1) and (2) the photon-photon correlation function G

(2)

Class (r 1 , r 2 ) to G (2) Class (r 1 , r 2 ) = n2 A + 2 nA nB + n2 B + +2 nA nB cos (ϕ 2 -ϕ 1 ) , (3)
where the relative phase ϕ i is given by ϕ i ≡ ϕ(r i ) = k d sin(ξ(r i )) (cf. Fig. 1) [START_REF]Note that the initial phases ψA and ψB[END_REF]. The visibility

V := G (2) max -G (2) min G (2) max + G (2) min (4)
of this classical signal calculates to

V Class = 2 nA nB (n A + nB ) 2 .
(

) 5
If one beam has a much higher mean photon number than the other, e.g., nA nB , the visibility V Class goes to zero as 2 nB nA . On the other hand, for nA = nB , we obtain the maximal value V Class = 50% [5].

Using Eq. ( 1) for mixed quantum-classical or mixed quantum-semiquantum sources the photon-photon correlation function is found to be

G (2) (r 1 , r 2 ) = ( Â(-) 1 + B(-) 1 )( Â(-) 2 + B(-) 2 ) ( Â(+) 2 + B(+) 2 )( Â(+) 1 + B(+) 1 ) = Â(-) 1 Â(+) 1 B(-) 2 B(+) 2 + + B(-) 1 B(+) 1 Â(-) 2 Â(+) 2 + Â(-) 1 Â(+) 2 B(-) 2 B(+) 1 + B(-) 1 B(+) 2 Â(-) 2 Â(+) 1 + B(-) 1 B(-) 2 B(+) 2 B(+) 1 (6)
where we have used the fact that the two photon sources are uncorrelated. In difference to the pure classical case above we can omit in Eq. ( 6) the term

Â(-) 1 Â(-) 2 Â(+) 2 Â(+) 1
as it vanishes identically for single photon emitters. However, unlike the pure quantum mechanical signal [28] consisting of, e.g., two single photons scattered by two atoms, we additionally have to take into account in Eq. ( 6) the term B(-

) 1 B(-) 2 B(+) 2 B(+)
1 . This term represents the probability that source B has emitted two photons which are subsequently measured at the two detectors. In the following, we start to discuss the configuration where source B is a standard classical source (Sect. II A). Thereafter, we study the configuration where source B is either in a PACS or PATS, i.e., in a state at the quantum-classical border (Sect. (II B)). These states are particularly interesting as they display properties ranging from completely classical to completely quantum.

A. Source B as a classical photon source

coherent source

Let us start to consider source B to be in a coherent state |α . Expanded in terms of Fock states |n , these states take the well known form [11]

|α = e -|α| 2 /2 ∞ n=0 α n √ n! |n , (7)
where α is any complex number. The k-th moment of the photon number operator a † a in the coherent beam can be calculated to

a †k a k C = α|a †k a k |α = |α| 2 k = nk , (8)
where n = a † a ∝ I B ∝ |E B | 2 is the mean photon number with E B , the amplitude of the electric field of source B. Using Eq. ( 6), the photon-photon correlation function of one quantum and one coherent source is then found to be

G (2) C (ϕ 1 , ϕ 2 ) = n2 + 2 n (1 + cos (ϕ 2 -ϕ 1 )) , (9)
so that the visibility of the mixed quantum-coherent photon-photon correlation function G

(2)

C (ϕ 1 , ϕ 2 ) becomes V C = 1 1 + n 2 . ( 10
)
For a classical coherent plane wave with n 1 we obtain V C = 100%, for n = 1 we arrive at V C = 66%, and for n → ∞ we have V C = 0. In [24] a single photon created by SPDC was entangled with an attenuated laser beam on a beam splitter. The authors obtained a modulation depth of (84 ± 3.2)%. According to Eq. (10) this corresponds to a mean photon number in the coherent beam of 0.293 < n < 0.475.

For the moments a †k a k T one calculates

a †k a k T = T r a †k a k ρT = k! nk , (12)
With these expressions Eq. ( 6) takes the form

G (2) T (ϕ 1 , ϕ 2 ) = 2 n2 + 2 n (1 + cos (ϕ 2 -ϕ 1 )) . ( 13
)
The visibility of the mixed quantum-thermal signal G

T (ϕ 1 , ϕ 2 ) thus calculates to

V T = 1 1 + n , (14)
what leads again to V T = 100% for n 1 and to V T = 0 for n → ∞. If we choose n = 1 we obtain V T = 50%, corresponding to the visibility of a pure classical photon-photon correlation signal exhibited by coherent light fields if nA = nB (cf. Eq. ( 5)). In a recent experiment a visibility of 82% was measured in this configuration, equivalent to a mean photon number in the thermal beam of n = 0.22 [26].

B. Source B as a semiquantum photon source

PACS

Now, let us consider the case where the field of source B is in a photon added coherent state [13]. In terms of the coherent states |α a normalized single photon added coherent state (in the following abbreviated P C) can be written as

|α, 1 = a † |α √ 1 + n . (15)
Hereby, n corresponds to the mean photon number in the coherent part of the light field. The k-th moments of the photon number operator a † a then take the form

α, 1|a †k a k |α, 1 = α|a a †k a k a † |α 1 + n , (16)
so that the first and the second moments read

n P C ≡ a † a P C = n2 + 3 n + 1 1 + n , a †2 a 2 P C = n2 + 4 n . ( 17
)
where we have introduced the net photon number n P C of the PACS field. With the help of Eq. ( 6), we arrive at

G (2) P C (ϕ 1 , ϕ 2 ) = n2 + 4 n + + 2 n2 + 3 n + 1 1 + n (1 + cos (ϕ 2 -ϕ 1 )) , (18)
and the visibility V P C thus takes the form

V P C = 1 1 + n3 +5 n2 +4 n 2 (n 2 +3 n+1) . ( 19
)
In case of a mixed quantum-coherent or mixed quantum-thermal photon-photon correlation signal (cf. Eqs. ( 9) and ( 13)) the mean photon number n in the coherent or thermal beam of source B trivially corresponds to the net photon number of the field, n C = n and n T = n, respectively. However, if we want to compare the visibility V P C of the mixed quantum-PACS photon-photon correlation signal to the foregoing results (cf. Eqs. ( 10) and ( 14)), we have to express V P C in terms of the net photon number n P C . We can rewrite Eq. ( 19) by using the identity (cf. Eq. ( 17))

n = - 3 2 + 1 2 n P C + + 1 2 5 -2n P C + n 2 P C (20)
and arrive at a visibility V P C depending on n P C rather than n. Since the analytic expression is rather complex, we do not present the explicit result -however, the corresponding outcome is plotted in Fig. 2.

For n P C → 1 we obtain V P C = 100%, whereas for n P C → ∞ we have V P C = 0, like in the mixed quantumclassical cases above. However, the visibility V P C of the mixed quantum-PACS photon-photon correlation function always displays higher values for a given net photon number a † a than the visibility V C for a mixed quantum-coherent source. For example, for a net photon number n P C = 1, we obtain a visibility V P C = 100% in contrast to a visibility V C = 66% for a net photon number n C = 1 in case of a mixed quantum-coherent source.

PATS

Next, we consider the field of source B to be in a PATS [14]. In case of a single photon added thermal state (in the following abbreviated as ST), the normalized density operator ρPT can be written in the Fock basis as [START_REF] Bell | Speakable and Unspeakable in Quantum Mechanics[END_REF] where n corresponds to the mean photon number of the thermal part of the light field. In contrast to the thermal density operator of Eq. ( 11) it is obvious that the vacuum term is missing and higher excited terms are rescaled. For the first and second moments of the photon number operator a † a we obtain

ρPT = 1 n + 1 a † ρT a = = 1 n (n + 1) ∞ n=0 n n + 1 n n |n n| ,
n P T ≡ a † a P T = 2 n + 1 , a †2 a 2 P T = 6 n2 + 4 n ( 22
)
where we have introduced again the net photon number n P T of the PATS field of source B. Thus Eq. ( 6) calculates to

G

(2)

P T (ϕ 1 , ϕ 2 ) = 6 n3 + 10 n2 + 4 n n + 1 + ( 2 3 ) + 2 (2 n + 1)(1 + cos (ϕ 2 -ϕ 1 )) ,
and for the visibility V P T we derive

V P T = 1 1 + 6 n3 +10 n2 +4 n 2 (n+1)(2 n+1) . (24)
Again, we have to formulate the visibility V P T in terms of the net photon number n P T in order to compare it to the previous mixed quantum-classical results. Inverting Eq. ( 22)

n = n P T -1 2 (25)
we obtain

V P T = 4 n P T 3 n 2 P T + 2 n P T -1 . (26)
From Eq. ( 26) we can see that for n P T → 1 the visibility V P T goes to 100% whereas for n P T → ∞ we have V P T = 0. Moreover, again the visibility V P T is taking higher values for any net photon number than V T of the mixed quantum-thermal state (cf. Eq. ( 14)). If we choose for example n P T = 2 we obtain V P T = 53.3%, whereas for n T = 2 we get V T = 1 3 (cf. Eqs. ( 12) and ( 14)). Fig. 2 displays the visibility of the various photonphoton correlation functions as a function of the net photon numbers n j (j = {T, P T, C, P C}) for the investigated combinations of mixed quantum-classical and mixed quantum-semiquantum sources. It thus shows a summary of the results discussed in this section. For comparison, the constant visibility V = 1 obtained for a pure quantum signal, e.g., two photons from two single photon emitters, and the constant visibility obtained for a pure classical signal (V Class = 1 2 ) are also plotted.

III. TEST OF THE QUANTUM CHARACTER OF THE MIXED RADIATION FIELDS

After deriving the various photon-photon correlation signals of the mixed quantum-classical and quantumsemiquantum sources in section II, the natural question arises of how to characterize and specify the nonclassicality [29,30] of these sources. For this purpose we introduce a particular version of the CS inequality for our photon-photon correlations, valid for classical fields. If this inequality is violated the underlying radiation field is non-classical. This is because the key ingredient in the derivation of the CS inequality is the assumption that the Glauber-Sudarshan P-function [11,31] behaves like a classical probability distribution. In its simplest form it reads [19]

a † 2 a 2 b † 2 b 2 ≥ a † b † a b 2 , ( 27
)
where a and b are the mode variables for the classical fields. For quantum fields, a and b are the annihilation operators. In this case the inequality can be violated since the Glauber-Sudarshan P-function can be negative, singular or need not exist. However, for the two point photon-photon correlations of our sources we have to derive an alternate form of the CS inequality. As discussed in Appendix A a direct application of Eq. ( 27) does not work.

Let us consider the net intensity I(ϕ) of two sources A and B at the point ϕ(r). If the underlying probability distribution (represented by the P-function) is positive, then the expression

: [α( Î(ϕ 1 ) -Î(ϕ 1 ) ) + β( Î(ϕ 2 ) -Î(ϕ 2 ) )] 2 : ( 28
)
should be positive for arbitrary α and β (where : : abbreviates normal ordering). As shown in the appendix the corresponding photon-photon correlations must satisfy this inequality provided that the P-function is positive. Any violation of this inequality implies then that the underlying radiation field is a nonclassical field. From Eq. ( 28) we arrive at the inequality (see Appendix A)

S := G(2) (ϕ 1 , ϕ 2 ) 2 G(2) (ϕ 1 , ϕ 1 ) G(2) (ϕ 2 , ϕ 2 ) ≤ 1 , ( 29
)
where the variance G(2) (ϕ k , ϕ l ) is given by

G(2) (ϕ k , ϕ l ) := G (2) (ϕ k , ϕ l ) -I A + I B 2 , (30) and G
(2) j (ϕ k , ϕ l ) (j = {T, P T, C, P C} corresponds to the photon-photon correlation functions given in Eqs. (3), ( 9), ( 13), ( 18) and (23). Hereby, I A (I B ) abbreviates the intensity of the photon sources located at R A (R B ) (cf. Fig. 1). Keeping in mind that we suppose uncorrelated sources and that we have I A QM = 1, and taking into account the corresponding moments given in Eqs. ( 8), ( 12), ( 16) and ( 22), the variances G(2)

j (ϕ k , ϕ l ) are calculated to G(2) Class (ϕ k , ϕ l ) = 2 n2 c kl , (31)
G(2) C (ϕ k , ϕ l ) = 2 n c kl -1 , (32)
G(2) T (ϕ k , ϕ l ) = n2 + 2 n c kl -1 , (33)
G(2) P C (ϕ k , ϕ l ) = -(3 n2 + 4 n + 2) (n + 1) 2 + + (2n 3 + 8n 2 + 8n + 2)c kl (n + 1) 2 , (34)
G(2) P T (ϕ k , ϕ l ) = 2 n2 + 2 (2n + 1)c kl -2 , (35)
G(2) (ϕ k , ϕ l ) = 2 c kl -2 , (36)
with the abbreviation c kl ≡ cos (ϕ k -ϕ l ). Note that for comparison we have also included the pure quantum mechanical case, i.e., two single photon emitters, where 29) and text).

G (2) (ϕ 1 , ϕ 2 ) = 2 (1 + cos (ϕ 2 -ϕ 1 )) (37)
with V = 1 [6]. We emphasize again that in case of a mixed quantum-PACS field or a mixed quantum-PATS field, we have to write G(2) P C and G(2) P T as a function of the net photon numbers n P C and n P T , respectively, to appropriately compare these signals to the mixed quantumclassical photon-photon correlation signals. To this end we have to insert the identities for the net photon numbers of Eq. ( 20) and (25) in Eqs. (34) and (35), respectively. If we now plug the variances of Eqs. ( 31)- (36) in the CS inequality Eq. ( 29), we obtain a violation for every mixed quantum-classical and mixed quantumsemiquantum signal below a certain net photon number a † a j (j = {T, P T, C, P C}), after an optimization with respect to ϕ 1 and ϕ 2 , i.e., with respect to the detector positions. Note that according to Fig. 2, reducing the net photon number n j in the field of source B increases the visibility V j of the corresponding photon-photon correlation signal (j = {T, P T, C, P C}).

The corresponding behavior of the Schwarz function S max j as a function of the visibilities V j (j = {T, P T, C, P C}) is displayed in Fig. 3. With the pure classical signal it is not possible to violate the CS inequality as S max Class = 1 for all V Class . For a mixed quantumthermal signal (T) the CS inequality is only violated for a visibility V T > 50%, i.e., if the net photon number in the thermal beam fulfills n < 1 (cf. Eq. 14). Note that a visibility > 50% is just above the maximal visibility of the pure classical coherent signal (cf. Eq. ( 5)). Moreover, for V T = 1/ √ 2 the CS function S max T diverges. By contrast, if we consider the mixed quantum-PATS signal (PT), the CS inequality can be violated already for a visibility V P T > 37.5%, corresponding to a net photon number n P T < 3 of source B (cf. Eq. ( 26)). Furthermore, we have S max P T → ∞ for V P T → 1. Considering the mixed quantum-coherent signal (C) and the mixed quantum-PACS signal (PC), if we choose a phase difference of π between ϕ 1 and ϕ 2 , we have a violation of the CS inequality for any finite net photon number of source B, what corresponds to a required visibility of these signals > 0 (see Fig. 2 and Eqs. (10), (19) and ( 20)). For the mixed quantum-coherent signal (C) we find S max C → ∞ for a visibility V C = 80%, corresponding to a net photon number of n = 1/2 in the coherent beam (cf. Eq. ( 10)). For the mixed quantum-PACS signal (PC) the CS function goes to ∞ for V P T → 1 (see Fig. 2). In this way we can pin down quantitatively the transgression to the non-classical regime for all our investigated mixed quantum-classical and mixed quantum-semiquantum signals. As expected, a signal of two single photon sources also violates the CS inequality, what is not depicted in Fig. 3. Finally it should be borne in mind that from no violation of the CS inequality no conclusion can be drawn regarding the quantum nature of the radiation field.

IV. TEST OF THE NON-LOCAL CHARACTER OF THE MIXED RADIATION FIELDS

In the past there have been many tests of Bell inequalities using photons obtained from parametric down converters (e.g., [32][33][34][35][36]). Further Bell inequalities have been studied using photons obtained either from two independent down converters [34] or one photon from a downconverter and one photon from a coherent source [24]. In all these experiments Bell inequalities were tested using polarization (spin like) correlations between the photons. We have already seen in the previous section that the radiation field produced by two disparate sources, where one source is a single photon emitter, is strictly quantum in nature in certain ranges of the net photon number of source B. Hence we would expect a violation of Bell-type inequalities [37]. However, in order to see this in detail we need to normalize our photon-photon correlation funcion. We use the following normalization (see Appendix B)

g(2) (ϕ 1 , ϕ 2 ) = G (2) (ϕ 1 , ϕ 2 ) N -1 , (38)
where N is an appropriate scaling factor. This scaling factor can be obtained by physical considerations -we take it to be the integrated photon-photon correlation function obtained by fixing one detector and moving the other detector. Using the moments given in Eqs. ( 8), ( 12), ( 16) and [START_REF] Clauser | [END_REF], and recalling that we assume uncorrelated sources and that we have I 2 A QM = 0, we can express the normalized photon-photon correlation functions in the collective form

g(2) j (ϕ 1 , ϕ 2 ) = G (2) j (ϕ 1 , ϕ 2 ) 2 π -1 G (2) j (ϕ 1 , ϕ 2 )dϕ 1 | ϕ2=const. -1 = = (I A + I B ) 2 + 2 I A I B cos (ϕ 2 -ϕ 1 ) (I A + I B ) 2 -1 = = V j cos (ϕ 2 -ϕ 1 ) , (39)
with the visibilities V j given in Eqs. ( 5), ( 10), ( 14), ( 19) and ( 26), respectively. We want to emphasize that our proposed rescaled correlation functions are independent of experimental inefficiencies, e.g., the quantum efficiency of the detectors and the restricted angle subtended by the detector surfaces (see Appendix B). Eq. ( 39) is the celebrated correlation commonly used in testing violations of Bell-type inequalities [START_REF] Bell | Foundations of Quantum Mechanics[END_REF]. This is also the correlation for two spins in a Werner state which is a mixed entangled state [41]. It is well known that such a correlation violates Bell's inequalities in the CHSH formulation [START_REF] Clauser | [END_REF] |g

j (ϕ 1 , ϕ 2 ) - g(2) j (ϕ 1 , ϕ 2 ) + + g(2) j (ϕ 1 , ϕ 2 ) + g(2) j (ϕ 1 , ϕ 2 )| ≤ 2 (40) if V j > 1 √ 2 . (2)
Down to the present day this particular version of Bell's inequalities is tested in most experiments (cf., e.g., [24,[33][34][35]). From Fig. 2 and Eqs. ( 10), ( 14), (19), (20) and (26), we can derive the net photon numbers n j of source B for which the CHSH inequalities are violated: The net photon number must be less than (a) √ 2 -1 for a thermal source, (b) ≈ 1.45 for a photon added thermal source (n ≈ 0.22), (c) 2 ( √ 2 -1) for a coherent source and (d) ≈ 1.52 for a photon added coherent source (n ≈ 0.29). Clearly in order to see the violations of locality the condition on n j for photon added sources is more relaxed.

V. SUMMARY

In conclusion we studied position dependent photonphoton correlations of photons emitted from disparate mixed quantum-classical and mixed quantumsemiquantum sources in free space. As quantum source we considered a spontaneously emitting atom; however we want to emphasize that our results are not limited to this case as any single photon source can act as a quantum source in the presented scheme [43]. The classical sources were represented either by thermal or coherent beams. These two classical sources have been studied in the past by overlapping photons at a beam splitter [24,26]. As semiquantum source we investigated both, PATS and PACS, what corresponds to a whole new class of disparate mixed quantum photon sources at the quantum-classical boundary. We introduced the CS inequality and obtained quantitative results for the nonclassicality of all our mixed radiation fields, depending on the visibility V j of the corresponding photon-photon correlations (see Fig. 3), what relates also to the net photon number n j in the field of source B (see Fig. 2). By normalizing the photon-photon correlation function in a physical way we obtained a violation of Bell-type inequalities for all considered mixed fields, again dependent on the visibility V j and the net photon number n j of source B. The violation of Bell-type inequalities requires more stringent conditions on the visibility of the photon-photon correlations than the violation of the CS inequality. This is reminiscent of the behavior of Werner states for spins.

  2. thermal sourceNow, let us assume that source B exhibits thermal statistics. The corresponding density operator expanded in Fock-states |n takes the form ρT = ∞ n=0 nn (n + 1) n+1 |n n| .

FIG. 2 :

 2 FIG. 2: Plot of the visibility of the photon-photon correlation function for different light sources dependent on the net photon number nj (j = {T, P T, C, P C}). QM (Class) abbreviates the visibility of a pure quantum (classical) source. Bell depicts the border to violate locality (cf. section IV for details).

FIG. 3 :

 3 FIG. 3: Plot of the maximum of the Schwarz function S max j dependent on the visibility Vj (j = {T, P T, C, P C}) for combinations of quantum and classical photon sources (cf. Eq. (29) and text).

  VI. ACKOWLEDGEMENTSR.W. gratefully acknowledges financial support by theElite Network of Bavaria. This work was supported by the Deutsche Forschungsgemeinschaft.

APPENDIX A: CS INEQUALITY FOR PHOTON-PHOTON CORRELATIONS

Let us investigate Eq. (28) in the form T r ρ [α δ Î(ϕ 1 ) + β δ Î(ϕ 2 )] 2 ≥ 0 ∀ α, β , (A1) with δ Î(ϕ i ) ≡ Î(ϕ i ) -Î(ϕ i ) and the density operator ρ in the Glauber-Sudarshan P-representation ρ = P (α)|α α|dRe(α)dIm(α) . (A2) P (α) abbreviates a (quasi-)probability function depending on the state of the investigated light field. For classical states of light P (α) is a classical probability distribution assuming only positive values. Since for any α and β the quantity [αδ Î(ϕ 1 ) + βδ Î(ϕ 2 )] 2 is also positive, inequality (A1) holds for any classical light field. However, for nonclassical states P (α) is a quasi-probability distribution and the P-function can be negative, singular or need not exist. Thus, the inequality can be violated. From Eq. (A1) it follows

Evaluating the determinant we arrive at

If we now use the identities

we obtain Eq. ( 29).

We want to emphasize that if instead of Eq. (A1) we use arguments similar to those leading to Eq. ( 27) we would get a CS inequality for normalized photon-photon correlation functions g (2) (ϕ i , ϕ j ) of the form (cf., e.g., [38])

For two quantum sources this form of the CS inequality is never violated as

Thus, clearly such an inequality does not enable us to test the quantum nature of our sources.

APPENDIX B: DERIVATION OF POSITION DEPENDENT BELL CORRELATIONS

In this appendix we want to derive Eqs. (38) and [START_REF] Bell | Foundations of Quantum Mechanics[END_REF]. Bell inequalities in the CHSH formulation are inequalities based on a double channel experiment. In order to establish this double channel character, we have to slightly change our setup. Let us consider that the photon sources A and B are connected to both detectors via optical fibers. Each detector shall now consist of two photomultipliers, located at the output ports of a 50-50 beamsplitter, representing a transmitting + and a reflectingchannel. Measuring the photon-photon correlation function, we thus consider four photodetectors, say D + 1 , D - 1 , D + 2 and D - 2 . In this configuration we can make use of the well known correlation function [START_REF] Bell | Foundations of Quantum Mechanics[END_REF][START_REF] Grangier | [END_REF] commonly employed in Bell experiments

where α, β = {+, -}. Hereby, G

α β (ϕ 1 , ϕ 2 ) abbreviates the photon-photon correlation of measuring one photon in the α channel at position ϕ 1 (i.e., at detector D α 1 ) and the other photon in the β channel at position ϕ 2 (i.e., at detector D β 2 ). In an experiment these correlations would correspond to four twofold coincidence rates. Keeping in mind that due to reflection we get a phase shift of π we arrive at

) where the plus (minus) sign holds for α = β (α = β). Thus, C(ϕ 1 , ϕ 2 ) calculates to

with the visibility

Eq. (B3) is equivalent to Eq. ( 39), i.e., Eq. ( 39) integrates the foregoing procedure of deriving correlations suitable for a Bell test. We want to emphasize that this expression holds for every position dependent mixed photonphoton correlation function with visibility V j . Due to the normalization on α,β G

α β (ϕ 1 , ϕ 2 ) the correlation C(ϕ 1 , ϕ 2 ), like Eq. (38), becomes independent of experimental inefficiencies of the form η = α 0 ∆Ω 4π which we have set to unity in the course of this paper. Hereby, α 0 stands for the quantum efficiency of the detectors and ∆Ω for the restricted solid angle subtended by the detector surfaces.