Two-axis two-spin squeezed states

The states generated by the two-spin generalization of the two-axis countertwisting Hamiltonian are examined. We analyze the behavior at both short and long timescales, by calculating various quantities such as squeezing, spin expectation values, probability distributions, entanglement, Wigner functions, and Bell correlations. In the limit of large spin ensembles and short interaction times, the state can be described by a two-mode squeezed vacuum state; for qubits, Bell state entanglement is produced. We ﬁnd that the Hamiltonian approximately produces two types of spin Einstein-Podolsky-Rosen (EPR) states, and the time evolution produces aperiodic oscillations between them. In a similar way to the basis invariance of Bell states and two-mode squeezed vacuum states, the Fock state correlations of spin EPR states are basis invariant. We ﬁnd that it is possible to violate a Bell inequality with such states, although the violation diminishes with increasing ensemble size. Effective methods to detect entanglement are also proposed, and formulas for the optimal times to enhance various properties are calculated.


I. INTRODUCTION
The concept of quantum squeezing has been central to the development of quantum metrology and its applications in quantum information science [1][2][3][4].In a squeezed state, it is possible to reduce the quantum mechanical noise (as measured by the variance) of an operator at the expense of another operator, according to the Heisenberg uncertainty relation.The most common experimental realizations of squeezing remain in optical systems [5][6][7].In addition to squeezed states for one mode, the two-mode squeezed state [8] is one of the central ingredients of optical-based quantum technologies, due to the entanglement that is possessed by this state.In a two-mode squeezed state, a linear combination of variables involving the two modes experiences suppression of quantum noise [2,3,9,10].Observables in the two modes obey Einstein-Podolsky-Rosen (EPR) type correlations, which can be employed for various quantum information tasks [9].Squeezed states of light have found applications in gravitational wave detection [11], quantum computation [9], interferometry [12], quantum metrology [13], and quantum cryptography [14], to name a few examples.
Atomic gases are another system that exhibits quantum mechanical squeezing.This includes atomic ensembles and Bose-Einstein condensates (BECs), where the relevant degrees of freedom are the internal spin states of the atoms [15].In contrast to the quantum optical case where the system is described by bosonic modes, for atomic gases the appropriate description is in terms of the collective spin of all the atoms.The seminal theoretical works of Kitagawa and Ueda [16] studied two main types of squeezed states in such atomic ensembles, the one-and two-axis spin squeezed states.These produce squeezing in the collective spin variables in a similar way to optical states.Such squeezed spin states have been observed in Refs.[17][18][19][20].
While squeezing on one spin has been studied in great detail theoretically and experimentally, the analog of twomode squeezing for the spin case is relatively less developed.The most widely known results for two spins are for atomic ensembles pioneered by the group of Polzik and co-workers [21][22][23].In these works, while the physical system is an atomic ensemble, the regime that is examined is where spin variables can be approximated by bosonic modes, according to the Holstein-Primakoff transformation.In these works the entangled state that is produced can be described within this approximation as a two-mode squeezed state.However, due to the fundamentally different nature of the spins to modes, it is well known that different dynamics can result for long evolution times, observed in effects such as oversqueezing [24][25][26][27][28][29].Entanglement between different spatial regions of a single BEC was experimentally observed [30][31][32].The two-spin version of one-axis squeezing was studied in Refs.[33,34] where it was found that a complex structure of entanglement is present, with a time dependence showing fractal characteristics.Methods to generate this state have been examined in numerous works [33][34][35][36][37][38][39].Procedures to generate other types of entangled states have also been studied [40,41].Such entangled states have been studied for various applications such as quantum computing [42,43] and quantum information [17,19,20,[44][45][46].
In this paper, we study the two-spin version of the twoaxis spin squeezed state, which we call the two-axis two-spin squeezed (2A2S) state.Within the Holstein-Primakoff approximation, corresponding to small interaction times, this state is equivalent to the two-mode squeezed state.In this approximation, the state has been studied before [47].Beyond these times, the dynamics starts to differ significantly from the two-mode squeezing interaction.We will particularly focus on understanding the states for longer interaction times, where the Holstein-Primakoff approximation is no longer valid.For the one-spin case, longer interaction times were studied in Refs.[26,48].We will generally be interested in the regime where the spin is large but finite, which is applicable for atomic gases.Even for BECs where the atomic numbers are much less than in thermal atomic ensembles, there can be easily more than 10 3 atoms in an ensemble making the spins very large [15,18,46].We study the two-spin two-axis spin squeezed state through various quantities, calculating correlations and probability distributions (Sec.III), entanglement (Sec.IV), Wigner functions (Sec.V), and Bell correlations (Sec.VI).In addition to elucidating the nature of the dynamics and the state, we interestingly find that the state can violate a Bell inequality without parity measurements.

II. TWO-AXIS TWO-SPIN SQUEEZED STATE
The system that we shall consider consists of two neutral atomic ensembles or BECs, where each atom has two relevant internal states.A common choice for the internal states are hyperfine ground states, such as the F = 1, m F = −1 and F = 2, m F = 1 states in the case of 87 Rb [15].In the case of BECs we denote the bosonic annihilation operator for the two states as a j , b j , respectively, where j ∈ {1, 2} labels the two BECs.These operators can be used to define an effective spin using the Schwinger boson operators where jkl is the completely antisymmetric tensor.For atomic ensembles, the total spin operators are written where σ k j,l is a Pauli operator for the lth atom in the jth ensemble.For simplicity, we take the total number of atoms N to be equal in both ensembles.The total spin operators also have the same commutation relation (2).For an atomic ensemble where the initial state and the Hamiltonian are completely symmetric with respect to particle interchange on a single ensemble, exactly the same results are obtained with either (1) or (3).Mathematically, the BEC form (1) is more convenient for calculations, and hence we will use this throughout this paper.However, it should be understood that our results equally apply to the atomic ensemble case.
The two-axis two-spin (2A2S) Hamiltonian is then defined as where and J is an energy constant.This is a straightforward generalization of the two-axis one-spin (2A1S) countertwisting Hamiltonian studied by Kitagawa and Ueda [16]: The 2A1S Hamiltonian produces squeezing and antisqueezing in the transformed variables respectively [16].We note that a similar generalization was studied previously to generalize the one-axis one-spin (1A1S) twisting Hamiltonian [49] to the one-axis two-spin (1A2S) Hamiltonian [33,34] The benefits of the 2A1S Hamiltonian are that it can attain a higher level of squeezing than the 1A1S Hamiltonian, and that the axes for optimal squeezing are fixed [16].The 2A2S squeezed states are produced by a unitary evolution according to the Hamiltonian (4) for a time t, where we have defined a dimensionless time τ = Jt/ h.The initial states are maximally polarized states in the S z direction, the same as for the 2A1S Hamiltonian.Here we defined the spin coherent states as where θ, φ are the angles on the Bloch sphere, and |vac is the vacuum containing no atoms, and j = 1, 2 label the BECs.We also define the Fock states as The Fock states are eigenstates of the S z operator according to For small times we may expand the exponential in Eq. ( 10) to second order, where the first terms in the expansion are where the states we have used are Fock states (12).All terms in the superposition have the same Fock numbers for the two BECs, hence, we immediately observe that the state is perfectly correlated in S z j .We can also see that for |τ | 1/N, the population of the a state will be small.A more precise criterion will be obtained in the next section.
In contrast to the 1A2S Hamiltonian, for which it is straightforward to find an analytical expression for the state at arbitrary evolution times, the 2A2S Hamiltonian cannot be diagonalized using a linear transformation of the bosonic operators.Only recently, static solutions for the 2A1S Hamiltonian were found [50].We thus resort to numerical methods to study the state and its properties.For pure state evolution, this involves evolving a vector of dimension (N + 1) 2 , hence, reasonably large systems can still be accessed with numerical methods.
In this paper we will focus on the theoretical state that is produced by the 2A2S Hamiltonian.Although we will not discuss the experimental procedures for generating the 2A2S Hamiltonian here, we note that several possibilities exist for producing the state, where techniques to produce the 2A1S Hamiltonian can be applied to the two-spin case.Several works have considered techniques to generate the 2A1S Hamiltonian [51][52][53].This can then be converted to the 2A2S Hamiltonian using a split-squeezing approach such as that given in Refs.[39,54].An alternative procedure is to generate the same correlations using optical interference methods [21,40].We plan to examine the experimental methods to produce the 2A2S Hamiltonian in a later work.

III. CORRELATIONS AND PROBABILITY DISTRIBUTIONS A. Holstein-Primakoff limit
To obtain some intuition about the state (10), let us first examine the state for small evolution times.For a system with fixed particle number N, the Holstein-Primakoff transformation [55] between spin operators to a single bosonic mode can be made according to In our model, the initial state of the system is the spin coherent state |0, 0 , in which all atoms in the two BECs occupy the respective internal state labeled by b.For N 1 and sufficiently small evolution times, the population of the a state satisfies a † a N.This allows us to approximate the Holstein-Primakoff-transformed spin operators as S z ≈ N.
Note that this approximation breaks down for larger evolution times due to the finite number of atoms N. As the Hamiltonian acts on the system, the occupation numbers of the a states increase while those of the b states deplete.At a certain time, which can be explicitly seen in the subsequent sections, the condition a † a N is no longer satisfied and the approximation fails.
Applying (16) to the 2A2S Hamiltonian, we have This is exactly the two-mode squeezing Hamiltonian [2,3,9] considered in quantum optics.The transformation of the mode operators is We can deduce the time for which the Holstein-Primakoff approximation is valid by evaluating the population of the a 1 and a 2 states.The population of the two states is always equal a † 1 a 1 = a † 2 a 2 and we obtain where in the last step we assumed Nτ 1. Demanding that a † 1 a 1 N, we have the criterion for the validity of the Holstein-Primakoff approximation as Now, let us define the canonical position and momentum operators as For the choice of phase between the two terms in Eq. ( 17), the relevant operators are those that are rotated by 45 • with respect to the quadrature axes where we used the definitions (7).The correlations for which the quantum noise is suppressed are then x1 + x2 and p1 − p2 [ 9,56].This can be seen by evaluating which become suppressed for large squeezing times.The corresponding antisqueezed variables are

B. EPR-type correlations
We now directly evaluate the correlations produced by the 2A2S Hamiltonian numerically, without applying the Holstein-Primakoff approximation.From (23) we expect that the variances of the observables become suppressed, for short times when the Holstein-Primakoff approximation holds.The observables in the perpendicular directions ( 24) are the antisqueezed variables.In Fig. 1(a), the variances of the observables (25) are plotted for short timescales τ ∼ 1/N.We see that the two variances have exactly the same time dependence and take a minimum at a time which we define as the optimal squeezing time τ (sq)  opt .In the time region 0 τ τ (sq) opt , the variance agrees well with Holstein-Primakoff approximation, giving which follows from the relations (23) and the fact that Var(O sq , τ = 0) = 2N.Beyond these times the variance increases and no longer follows (27).For longer timescales, as shown in Fig. 1(b), the variance follows aperiodic oscillations between low and high variance states.Some relatively low variance states are achieved (e.g., particularly around τ ≈ 3), although the minimum variance at the times τ (sq)  opt is not attained again.
The antisqueezed variables (26) are shown in Fig. 1(c) for short timescales τ ∼ 1/N.Again, the two variables (26) have exactly the same time dependence, and initially increase according to which follows from (24).In contrast to genuine two-mode squeezing, the variance does not increase unboundedly but reaches a maximum.We call this time the optimal antisqueezing time τ (asq) opt , which we find is not exactly at the same time as the optimal squeezing time τ (sq)  opt .As with the squeezed variables, for longer timescales, as shown in Fig. 1(d), the antisqueezed variables show aperiodic oscillations, with a similar range to Fig. 1(b).Some low variances states are also present with the antisqueezed variables, again around τ ≈ 3.However, these states do not reach the low level of the variance attained with the squeezed variables (25).
The existence of a bound on the amount of squeezing, as opposed to unbounded genuine two-mode squeezing, is a consequence of the finite atom number N. In the large-N limit, the Holstein-Primakoff approximation is valid and our system is equivalent to two-mode squeezing.Restricting the atom number N to a finite value renders the Fock space finite dimensional and leads to oscillations of the populations between the north and south poles of the Bloch sphere, as will be seen in the following section.

C. Expectation values
It is also instructive to examine the expectation values of the spin operators in the 2A2S squeezed state.Figure 2 shows the expectation values of the operators S x j , S y j , S z j .Due to the symmetry between the initial state of the two ensembles and the 2A2S Hamiltonian, identical values are obtained for the two ensembles j = 1, 2. Furthermore, the expectation values of two of the operators are always zero: This can be seen from ( 14), where the Hamiltonian creates pairs of equal number Fock states.Since the S x and S y  (31).The fit parameters for minimized squeezing are p 0 = 0.467, p 1 = 0.508, maximum antisqueezing are p 0 = 0.700, p 1 = 0.530, zeros of S z are p 0 = 0.727, p 1 = 0.536, the fidelities with the spin EPR state are p 0 = 0.803, p 1 = 0.544.operators shift the Fock states by one unit, The expectation values of S x and S y are zero for all time.Meanwhile, the expectation value of S z j undergoes aperiodic oscillations and flips sign numerous times during the evolution.In particular, a sign change is observed in the vicinity of τ (asq)  opt .This can be understood from (14), where at τ ∼ 1/N the sum contains all terms with a similar magnitude.Again, the time where S z j = 0 is not exactly the same as the optimal squeezing time τ (sq)  opt or antisqueezing time τ (asq) opt .We label the first time that a sign flip of S z j occurs by τ (S z ) opt .We may thus picture the type of state that is produced as a two-spin version of the planar squeezed state, where the squeezing is observed in the (S x j , S y j ) plane [57].

D. Optimal squeezing times
In order to maximize the correlations between the two spins, it would be useful to have a general expression that gives the optimal squeezing time τ (sq)  opt .Unfortunately, we have not been able to find an analytical expression to give this for general N.However, it is possible to find approximate formulas that can give this to good accuracy.We will investigate this in the following.
As can be seen from Figs. 1(a), 1(c), and 2(a), the optimal squeezing time τ (sq)  opt , optimal antisqueezing time τ (asq) opt , and the time τ (S z ) opt of the first zero of S z do not necessarily coincide.The optimal times according to each criterion are shown in Fig. 3(a).Generally, the optimal squeezing time tends to give slightly smaller values than that for antisqueezing, or the time for the first zero of S z .The latter two conditions give very similar values of the optimal time.In Fig. 1(a) we notice that the time that the Holstein-Primakoff approximation starts to deviate from the exact expression coincides with the optimal time for squeezing.This suggests that the functional form of (20) may be a suitable form to fit the optimal times in Fig. 3(a).Using the fitting form we fit the various optimal times as shown in Fig. 3(b).This fitting function works very well, with the logarithmic term accounting for the nonlinear behavior seen in the 1/N plot.
For the optimal time extracted for the squeezed variables we obtain fit parameters p 0 = 0.467, p 1 = 0.508, close to the theoretically calculated values of (20), p 0 = ln 2, p 1 = 1 2 .Since ( 31) is guaranteed to approach τ opt → 0 in the limit of large N, we expect that this fit will interpolate the larger N optimal times with good accuracy.The parameters for the other optimal times are given in the caption of Fig. 3.
Using the optimal squeezing times, the maximal squeezing that can be attained can be estimated from the Holstein-Primakoff relation according to where in the last line we substituted the expression (31).The minimal squeezing level tends to improve approximately as 1/N for larger ensemble sizes due to the factor of N 2p 1 in the denominator, relative to the spin coherent state variance 2N.

E. Probability density of the two-axis cosqueezed state
Another way to visualize the correlations is to plot the probability distributions when the state (10) is measured in various bases.Specifically, we consider the Fock states which are the eigenstates of the rotated operators (7): These Fock states can be transformed from the S z -eigenstate Fock states (12) using the relations in Appendix A. The probability of a measurement outcome k 1 , k 2 for various Fock states is then where l 1 , l 2 ∈ {x, y, z}.
The probabilities for two evolution times before and near τ (sq)  opt are shown in Fig. 4. The effect of the correlations is seen in the ( Sx  measurement, we see Fock state correlations arising from the fact that the 2A2S Hamiltonian always produces Fock states in pairs, as shown in Eq. (14).
Comparing the ( Sx 1 , Sx 2 ) and ( Sy 1 , Sy 2 ) in more detail reveals some interesting effects.For the suboptimal squeezing time in Fig. 4(a) we see that the peak of the distribution is at Sx 1 = Sx 2 = 0, which is the original position of the Gaussian before the Hamiltonian is applied.In Fig. 4(b), we see that the maximum is at the ends of the distribution.Thus, the optimal squeezing corresponds approximately to tuning the time such that the ( Sx 1 , Sx 2 ) and ( Sy 1 , Sy 2 ) distributions are at their flattest.opt − τ (sq,asq, S z ) opt , where τ (E)  opt is the optimal time according to von Neumann entropy, τ (sq)  opt is for the squeezed variables, τ (asq) opt is for the antisqueezed variables, τ (S z ) opt is for zero of S z .
For the remaining correlation pairs, the distributions are always symmetrical in the variables Sx and Sy , hence, give zero when averaged.Thus, there is no correlation between the remaining variables.The lack of correlations in the off-diagonal combinations in Fig. 4 is also a feature of standard Bell states.The type of state that is considered here is therefore a natural generalization of the EPR correlated state for spin ensembles, which we discuss further in Sec.IV B.

IV. ENTANGLEMENT A. Time dependence of entanglement for pure states
We now turn to the entanglement that is generated in 2A2S state.The entanglement that we consider is that present between the two BECs, which forms a natural bipartition in the system.We will not consider other types of intraensemble entanglement here, as have been considered for single ensemble squeezed states [15].For the 2A2S squeezed state, there is no squeezing on a single BEC, hence, we do not expect intraensemble entanglement to be a relevant quantity in this context.
For pure states, bipartite entanglement can be quantified using the von Neumann entropy where is the reduced density matrix for BEC 2. Figures 5(a We see that the entanglement first reaches a maximum at a similar time to the optimal squeezing time τ (sq) opt , and reaches nearly the maximum possible entanglement between the two BECs.For larger values of N, the oscillations have a higher frequency, with a period that is ∼2τ opt .Figure 5(c) shows the maximal entanglement as a function of N. For each N, we find the maximum value of the entanglement by optimizing the time in the vicinity of the first maximum.We call the time that this occurs the optimal entanglement time τ (E )  opt .We see that the optimized entanglement approaches the maximum possible entanglement E max for large N.This is in agreement with estimates from the Holstein-Primakoff approximated Hamiltonian, where the maximum value is reached for N → ∞.The convergence to the maximum value, however, does occur logarithmically (see Appendix B), as observed by the slow approach of Fig. 5(c).The entanglement oscillates between large and small values and tends to occur at the values corresponding to S z = 0.This is true not only in the vicinity of the first maximum that is reached at τ (S z ) opt , but for all τ .Figures 5(a) and 5(b) mark all the times (with a dot in the figure) where S z = 0. We see that each peak in the entanglement occurs when S z = 0.This is reasonable to expect from the point of view that zero values of S z correspond to states that potentially have large degrees of correlation in the S x and S y variables.In Fig. 5(d) we compare the various optimal squeezing times to the time to maximize the entanglement τ (E )  opt .The time that most closely approximates the maximal entanglement is τ (S z ) opt , which slightly overestimates the optimal entangling time, but gives the closest approximation.We point out that S z = 0 of course does not ensure that entanglement is present, particularly for mixed states, since this can equally occur for dephased states.However, this is a convenient heuristic that could be easily measured that coincides with large values of entanglement.

B. Spin EPR state
We have seen in Fig. 5(c) that near-maximal entanglement can be obtained at optimized evolution times of the 2A2S Hamiltonian.We have also seen in Fig. 4(b) that at the optimized squeezing times, very flat distributions of the correlations can be obtained.These facts suggest that a good approximation for the state in the large-N regime is where we defined the state This state has the maximum possible entanglement E max between the two BECs, and exhibits squeezing in the variable Sx 1 + Sx 2 .Algebraic manipulation allows one to rewrite this state equally as which have the Sy 1 − Sy 2 and Sz 1 − Sz 2 correlations, in agreement with Fig. 4.Such a state is a type of spin EPR state which exhibits correlations in a similar way to Bell states and continuous variable two-mode squeezed states, in all possible bases [9].In fact, as we show in Appendix C, the correlations are for any choice of basis such that where the Fock states rotated by a polar θ and azimuthal φ angles are defined by In Fig. 6 show the fidelity of the 2A2S squeezed state with reference to the spin EPR state, defined as We also plot the fidelity with respect to another spin EPR state defined without phases: where We see that the state attains high overlap with the |EPR − at a time τ (E ) opt as expected, and oscillates with peaks at similar times as the peaks in the antisqueezing parameters seen in Fig. 1(d).On comparison with Fig. 5(b), we see that every second peak in the entanglement corresponds to the peaks for the fidelity F − .The remaining peaks occur for the fidelity F + .The timing of the peaks in F + also matches the peaks in squeezing parameters in Fig. 1(b).This makes it clear that the effect of the 2A2S Hamiltonian is to first generate a state closely approximating |EPR − , which then subsequently evolves to |EPR + , and this cycle repeats itself in an aperiodic fashion.
In Fig. 6(b) we examine the scaling of the fidelity with N. We optimize the interaction time τ such as to maximize F − in the region of the optimal squeezing time.The time where the first maximum in F − is attained is defined as τ (F )  opt .The optimal time τ (F )  opt is again found to be most similar to τ (S z ) opt , but not precisely the same (see Fig. 3).The fidelity approaches F − ≈ 0.9 for the largest ensemble sizes that we examined (N = 160).While the data suggest a ∝1/N relation, we do not perform an extrapolation for large N since we cannot completely rule out that logarithmic corrections may exist due to the fidelity being a sensitive quantity for large systems.It does, however, appear that the dependence with N is rather weak and it is not unreasonable to expect that high fidelity with spin EPR states can be generated even for realistic BEC sizes.

C. Entanglement detection
The results of Fig. 5 clearly show that the 2A2S squeezed state is entangled for all times except for τ = 0.For pure states, the von Neumann entropy completely quantifies the amount of entanglement in a bipartite system.However, calculation of the von Neumann entropy relies on the availability of the complete wave function |ψ (t ) , which may be difficult to extract experimentally, particularly for large dimensional systems.Tomographic reconstruction of the full density matrix has a high overhead in terms of the number of measurements that need to be made.Another potential experimental constraint is that only particular types of measurements may be feasible.The most common types of measurement in the context of BECs are Fock state measurements.However, at present the atom-number resolution is a limitation, hence, quantities that are insensitive to single-atom fluctuations are preferred.Thus, the experimentally preferred quantities are low-order spin expectation values.In this section, we discuss criteria which only involve low-order spin expectation values and can show the presence of entanglement between the two BECs.
We compare three potential criteria that can be used to detect entanglement using low-order spin expectation values.The first is the Giovannetti-Mancini-Vitali-Tombesi (GMVT) criterion [58], which states that any separable state obeys Here, g x , g y are parameters to be optimized, and we have chosen operator combinations that attain small values near the optimal squeezing times τ (sq) opt .We find that the optimum values of the parameters are in this case g x = g y = 1.The second criterion is the Duan-Giedke-Cirac-Zoller (DGCZ) criterion [59], which for the type of correlations that are present in this state reads as which is true for any separable state.The third criterion is the Hofmann-Takeuchi (HT) criterion [60] based on local uncertainty relations for three spin operators which is true for any separable state.Violation of the inequalities ( 46)-( 48) signals the presence of entanglement.46), the Duan-Giedke-Cirac-Zoller (DGCZ) criterion (47), and the Hofmann-Takeuchi (HT) criterion (48).Entanglement between the two BECs is detected for all values below the separability bound, represented by the shaded area.
range.Initially, all three criteria detect entanglement successfully.However, the GMVT and DGCZ criteria fail at times near the optimal squeezing time τ (sq)  opt due to the fact that at these times S z j becomes small.The HT criterion does not fail in these regions since it does not rely on a comparison with S z j .For longer times [Fig.7(b)], all criteria are generally less successful at detecting entanglement, only detecting some specific time regions.For times which show strong squeezing, the entanglement criteria are able to detect entanglement.The most successful criterion is the HT criterion, hence, for this particular state (48) appears to be the method of choice for correlation-based entanglement detection.

A. Definitions
In this section, we further analyze the 2A2S squeezed state by visualizing its Wigner function.The Wigner function represents the state as a quasiprobability distribution on the Bloch sphere and for two-mode BECs can be defined as [61] where j = N/2, the spherical harmonics Y lm (θ, φ), and ρ lm given by Here, ρ is the density operator of the state being analyzed and j 1 m 1 ; j 2 m 2 |JM are the Clebsch-Gordan coefficients coupling two angular momentum eigenstates | jm .These states can be expressed in terms of the S z eigenstates using the relation The Wigner function as in Eq. ( 49) is defined for a two-mode BEC with fixed atom number and thus cannot be directly applied to our case of two BECs involving four modes.Instead, we calculate the marginal and conditional Wigner functions where tracing or projective operation is performed to obtain the Wigner functions on a single BEC.

B. Marginal Wigner functions
To calculate the marginal Wigner function, take the partial trace of the two-BEC state (10) over BEC 1 to obtain the reduced density matrix ρ 2 .The state (36) is then used to calculate the Wigner function according to (49). Figure 8 shows the marginal Wigner function for four different interaction times τ .The initial state [Fig.8(a)] starts as a Gaussian at the north pole of the Bloch sphere.As the 2A2S Hamiltonian is turned on [Fig.8(b)], the diameter of the Gaussian increases, until the distribution nearly covers the whole Bloch sphere at the optimal squeezing time τ (sq)  opt [Fig.8(c)].At this point the average value of all spin variables is small, in agreement with Fig. 2. For longer squeezing times, the probability distribution becomes more concentrated at the south pole of the Bloch sphere [Fig.8(d)], at which point S z turns negative as seen in Fig. 2. The flipping between the north and south poles continues for longer times, as seen in Fig. 2(b).The Wigner function is completely rotationally symmetric around the S z axis at all times.In a similar way to the two-mode squeezed state, there is no squeezing for a single BEC, and the squeezing only appears in variables involving both BECs.

C. Conditional Wigner functions
For the conditional Wigner function, we first project the two-BEC state (10) onto a Fock state (42) of one of the BECs (say, BEC 1), for a particular basis specified by (θ, φ).This would correspond to performing a Fock state measurement in a given basis on BEC 1, where a random collapse occurs, then examining the measured state of BEC 2. The resulting quantum state is where the projector onto a Fock state on BEC 1 is To calculate the conditional Wigner functions, the state ρ = |ψ k (θ, φ, τ ) ψ k (θ, φ, τ )| is substituted into (49).
The resulting conditional Wigner functions for projections in various bases are shown in Fig. 9. From the form of the spin EPR state (41), we expect that the resulting state is |ψ k (θ, φ, τ opt ) ≈ |k (θ,π−φ) 2 (54) on BEC 2. For projections with k = N [Figs.9(a) and 9(b)], the projected state appears as a Gaussian Wigner function, characteristic of a spin coherent state, centered at the angular parameters (θ, π − φ).This is in agreement with (54), using the fact that |k = N (θ,φ) = |θ, φ . ( Similarly, for a projection with k = 0 the resulting state is a spin coherent state centered at (π − θ, φ − π ), using the fact that In the case shown in Fig. 9(c), the resulting spin coherent state is more distorted than those shown in Figs.9(a) and 9(b) because parameters are chosen such that the final state is located at θ = 3π/4, φ = −π/2.From the (S z 1 , S z 2 ) plot in Fig. 4(b), we see that the probability distribution tends to diminish for negative S z , i.e., near the south pole of the Bloch sphere.We attribute the distortion of the distribution in Fig. 9(c) to the generation of an imperfect spin EPR state due to the relatively small ensemble sizes N = 10 considered here.
Finally, projecting on k = N − 1 [Fig.9(d)] produces a single-particle Fock state centered around (θ, π − φ).The resulting state shows a distribution that is similar to that of a single-particle Fock state, with a negative central region, surrounded by a positive ring.Again, deviations from the exact spin EPR state cause some differences to an ideal Fock state Wigner function.For larger ensembles, we expect that the distributions will more closely follow the distributions of the Fock states.

VI. BELL'S INEQUALITY
In this section, we show that the 2A2S squeezed state violates a Bell inequality using only experimentally accessible correlations of total spin operators (1).Recently, the violation of a multipartite Bell inequality has been experimentally achieved with a single two-mode BEC [62].The measurement of violations of a bipartite Bell inequality using two spatially separated BECs, however, has not been reported yet.Previous studies on 1A2S entangled states suggested that parity measurements (i.e., ±1 depending on whether |k is even or odd) are required for violating a Bell inequality [54].Parity measurements are currently experimentally challenging because they are sensitive at the single-atom level.Methods to violate a Bell inequality without the use of parity measurement are therefore of interest.
We use the Bell-CHSH inequality [63] for two observers with two local measurement choices.Every local realist theory must satisfy C = M (1)  1 M (1)   2 where 2 are local measurement choices on BEC 1 and BEC 2, respectively, and i, j ∈ {1, 2} label the two measurement choices.The measurement operators leading to a violation of (57) are given by Here the meaning of the sgn function is that the sign of the eigenvalues of the spin operators are taken: sgn( Sx cos θ + Sy sin θ ) = N k=0 sgn(2k − N )|k (θ,0) k| (θ,0) , (59) where |k (θ,φ) are the eigenstates of the rotated spin operator Sx 1 cos θ + Sy 1 sin θ and is defined in Eq. (42).Taking the sgn function of the rotated operators produces dichotomic (i.e., two-valued) measurement outcome as needed for the Clauser-Horne-Shimony-Holt (CHSH) inequality.But importantly, this way of mapping the collective spin measurement to a dichotomic observable does not require single-atom precision since the variation of the eigenvalues is split into two broad regions, k < N/2 and k > N/2.
The choices of the parameters θ (i) n and the interaction time τ should generally be optimized such that a maximal violation is obtained.First, we have found that the largest violations occur for interaction times τ ≈ τ (sq)  opt .This is reasonable from the point of view that the Bell-CHSH inequality (57) is correlation based, and the maximal squeezing tends to maximize the correlations.For the angular parameters, we have found that a strong violation can be found using the parameter choices  (62).Interaction times of τ = τ (sq)  opt are used throughout.
For the qubit case N = 1, the above basis choice reduces to the well-known optimal CHSH measurements by taking θ B = π/2.For larger values of N, the sign-binned operators (58) no longer coincide with the spin operators.We find that the Tsirelson bound can no longer be reached for N > 1, possibly due to the loss of information caused by the binning process.We do, however, find violations of (57) for angles θ B significantly smaller than π/2.Both the violations and the optimal angle θ B diminish with increasing N.
Figure 10(a) shows the optimized Bell-CHSH violations with respect to the angles θ B .The amount of violation is found to empirically have scaling that agrees very well with the relation The fact that the effect diminishes for large N is expected as our system can be described by continuous variable mode operators in the limit of large N, according to the Holstein-Primakoff approximation.For two-mode squeezed states, it is known that no Bell violation is possible with measurements that are linear combinations of quadrature operators [9].Maximal violations of the CHSH inequality with EPR states have been reported with parity measurements, e.g., in Ref. [64], suggesting that single-atom resolution measurements may be needed to achieve maximal violations with the spin EPR state.
In the limit of large N we therefore expect that the Bell-CHSH expression with the coarse-grained measurement operators (58) approaches 2, which can be attained if θ B = 0, such that all four correlators are sgn(S x 1 S x 2 ) → 1. Interpolating small and large N, and from the highly linear nature of Fig. 10(a), we expect that violations can be found for any finite value of N, although they will become smaller and increasingly challenging to measure for larger N.
The optimal angles to violate the Bell-CHSH inequalities are given in Fig. 10(b).The angle does not follow a simple power-law relation, and thus we fit it with a Padé approximant which gives a close approximation to the numerically obtained angle.
We have examined the 2A2S squeezed state from multiple points of view: squeezing, spin expectation values, probability distributions, entanglement, Wigner functions, and Bell correlations.A consistent picture emerges from studying various quantities.Starting from two BECs which are polarized in the positive S z direction, the 2A2S Hamiltonian produces the spin EPR state |EPR − , the analog of the two-mode squeezed state for spins.This state then evolves to the state of two BECs polarized in the negative S z direction, from which then the spin EPR state |EPR + is produced.After returning to two BECs polarized in the positive S z direction, the process repeats itself.The picture that we describe here is not exact, even for the decoherence-free case that was examined in this paper.As can be seen in Figs. 1, 2, 5, and 6, the oscillations are not perfectly periodic in amplitude or period.For such imperfect oscillations one may expect that the oscillations die out relatively quickly.However, the oscillations are remarkably persistent for the times that we have examined, and the fidelities to the spin EPR states remain high after a large number of oscillations.
The optimal times to produce the largest amounts of squeezing, entanglement, Bell correlations, and fidelities with the spin EPR states are found to be similar, but not precisely the same.Thus, each quantity must be optimized separately in order to obtain the optimal values.Approximate formulas for the optimal times were obtained using a suitable fitting function which should be accurate particularly for large N since it is determined by interpolating between numerically determined intermediate N and N → ∞.It is found that optimal times for the entanglement and fidelities with spin EPR states are closely approximated by the times that S z = 0. Meanwhile the optimal times to violate the Bell-CHSH inequality are closest to the optimal squeezing time.Since the measurement of S z is relatively simple, this is a convenient heuristic that can be used to generate the desired state.
One interesting feature of the 2A2S squeezed states is that they are able to violate a Bell-CHSH inequality, and it appears that this is possible for all finite N. We deduce this from the fact that in the limit of large N, the 2A2S squeezed state should approach a two-mode squeezed state.The level of violation is rather small, dropping off as 1/N, which makes it more difficult to observe for large ensembles.However, we conjecture that this may be the best that can be done for a two-measurement, two-outcome CHSH inequality that only uses first-order correlators of coarse-grained collective spin operators since it is known that quadrature operator measurements alone cannot detect nonlocality in two-mode squeezed states [9].The Fock states of the spin operators corresponding to Sx and Sy are given by |k (x) = e −iS z π/8 e −iS y π/4 |k (z) , |k (ỹ) = e −iS z 3π/8 e −iS y π/4 |k (z) . (A1) Here the matrix elements of the S y rotation are given by where |k = |k (z) .

APPENDIX B: ENTANGLEMENT SCALING UNDER THE HOLSTEIN-PRIMAKOFF APPROXIMATION
In this Appendix, we derive the entanglement at the optimal times under the Holstein-Primakoff approximation.Starting from the Holstein-Primakoff approximated Hamiltonian (17), the two-mode squeezed wave unction is given as [9] e −iH 2 t |0 = c where in the second line we further approximated the expression for large N. We see that for large N, E /E max logarithmically approaches 1.

FIG. 2 .
FIG. 2. Expectation values of spin operators for the two-axis two-spin squeezed state for (a) short timescales τ ∼ 1/N; (b) long timescales τ ∼ 1.The number of atoms per ensemble is taken as N = 20.

FIG. 5 .
FIG. 5. Entanglement of the two-axis two-spin squeezed state, as measured by the von Neumann entropy E normalized to the maximum value E max = log 2 (N + 1), as a function of the interaction time τ .The number of atoms in each BEC is (a) N = 10, (b) N = 20.The dots in each figure denote the times when S z 1 = 0 for each N. (c) The maximum value of the von Neumann entropy as a function of N. (d) Difference in the optimal time τ opt = τ (E) opt − τ (sq,asq, S z ) ) and 5(b) show the von Neumann entropy normalized to the maximum value E max = log 2 (N + 1) for two N + 1 level systems.

FIG. 6 .
FIG.6.(a) Fidelities of the two-axis two-spin squeezed state(10) with respect to the spin EPR states(40) and(45).Time dependence of the fidelity for N = 20.(b) Optimized fidelity(40) as a function of 1/N.Inset shows zoomed-in region for large N.

2 FIG. 10 .
FIG. 10.(a) Optimized violation of the Bell-CHSH inequality C-2 as a function of 1/N for the 2A2S squeezed state.Dots show the numerically calculated values and the solid line is the linear fit (61).(b) Optimal Bell angles θ B with the basis choices (60) as function of 1/N.Dots show the numerically calculated values and the solid line is the Padé approximant(62).Interaction times of τ = τ(sq)  opt are used throughout.

ACKNOWLEDGMENTS
This work is supported by the Shanghai Research Challenge Fund; New York University Global Seed Grants for Collaborative Research; National Natural Science Foundation of China (Grants No. 61571301 and No. D1210036A); the NSFC Research Fund for International Young Scientists (Grants No. 11650110425 and No. 11850410426); NYU-ECNU Institute of Physics at NYU Shanghai; the Science and Technology Commission of Shanghai Municipality (Grants No. 17ZR1443600 and No. 19XD1423000); the China Science and Technology Exchange Center (Grant No. NGA-16-001); and the NSFC-RFBR Collaborative grant (Grant No. 81811530112).APPENDIX A: EXPRESSION FOR PROBABILITY OF TWO-AXIS TWO-SPIN STATE MEASURED IN VARIOUS SPIN BASES