Light mediated non-Gaussian atomic ensemble entanglement

TL;DR

This paper extends the analysis of light-mediated entanglement between atomic ensembles by solving the exact time evolution without the Holstein-Primakoff approximation, revealing non-Gaussian states at longer interaction times.

Contribution

It provides an exact solution for the entanglement dynamics, showing the transition from Gaussian to non-Gaussian states over time, beyond the standard approximation.

Findings

Short entangling times produce two-mode squeezed states.

Longer times lead to non-Gaussian entangled states.

Exact evolution reveals deviations from Holstein-Primakoff predictions.

Abstract

We analyze a similar scheme for producing light-mediated entanglement between atomic ensembles, as first realized by Julsgaard, Kozhekin and Polzik [Nature {\bf 413}, 400 (2001)]. In the standard approach to modeling the scheme, a Holstein-Primakoff approximation is made, where the atomic ensembles are treated as bosonic modes, and is only valid for short interaction times. In this paper, we solve the time evolution without this approximation, which extends the region of validity of the interaction time. For short entangling times, we find this produces a state with similar characteristics as a two-mode squeezed state, in agreement with standard predictions. For long entangling times, the state evolves into a non-Gaussian form, and the two-mode squeezed state characteristics start to diminish. This is attributed to more exotic types of entangled states being generated. We characterize the states by examining the Fock state probability distributions, Husimi $Q$ distributions, and non-local entanglement between the ensembles. We compare and connect several quantities obtained using the Holstein-Primakoff approach and our exact time evolution methods.

Equations

$$x=\frac{J^y}{\sqrt{2N}},p=\frac{J^z}{\sqrt{2N}},$$

$$H=g{S}^z(J_1^z+J_2^z)$$

$$\hslash \frac{d{x}_{\text{ph}}}{dt}=ig[p_{\text{ph}}(p_{\text{at},1}+p_{\text{at},2}),x_{\text{ph}}]=g(p_{\text{at},1}+p_{\text{at},2}),$$

$$\hslash \frac{d{p}_{\text{ph}}}{dt}=ig[p_{\text{ph}}(p_{\text{at},1}+p_{\text{at},2}),p_{\text{ph}}]=0.$$

$$x_{\text{ph}}(t)$$

$$p_{\text{ph}}(t)$$

$$x_{\text{at},1}(t)$$

$$x_{\text{at},2}(t)$$

$$p_{\text{at},1}(t)$$

$$p_{\text{at},2}(t)$$

$$p_{\text{at},1}(t)+p_{\text{at},2}(t)=p_{\text{at},1}(0)+p_{\text{at},2}(0)$$

$$x_{\text{at},1}(t)+x_{\text{at},2}(t)=x_{\text{at},1}(0)+x_{\text{at},2}(0)$$

$$J^x$$

$$J^y$$

$$J^z$$

$$a$$

$$b$$

$$J^x$$

$$J^y$$

$$J^z$$

$$|J^x=-N⟩=\frac{1}{\sqrt{N!}}(b_x^{†}{)}^N|0⟩.$$

$$J^x$$

$$J^y$$

$$J^z$$

$$|\psi ⟩=\sum _k{\psi }_k|k⟩$$

$$|k⟩=\frac{(a^{†}{)}^k(b^{†}{)}^{(N-k)}}{\sqrt{k!(N-k)!}}|0⟩$$

$$|k_x⟩=\frac{(a_x^{†}{)}^k(b_x^{†}{)}^{(N-k)}}{\sqrt{k!(N-k)!}}|0⟩.$$

$$|\psi ⟩=\sum _{k_x}{\psi }_{k_x}^{\prime }|k_x⟩$$

$$|k_x⟩\text{ (spins)}\leftrightarrow |k_{\text{ph}}⟩\equiv \frac{1}{\sqrt{k!}}(a_x^{†}{)}^k|0⟩\text{ (mode)}.$$

$$|\psi_{\text{light}}(t=0)\rangle=|\xi\rangle=e^{-\frac{|\xi|^{2}}{2}}\exp\Big{[}\frac{\xi}{\sqrt{2}}(c^{\dagger}+d^{\dagger})\Big{]}|0\rangle$$

$$S^x$$

$$S^y$$

$$S^z$$

$$|{\psi }_{\text{atoms}}(t=0)⟩=|\frac{\pi }{2},0⟩{⟩}_1|\frac{\pi }{2},\pi ⟩{⟩}_2$$

$$|\theta ,\varphi ⟩{⟩}_j=\prod _{n=1}^{N_j}\left(\cos \left(\frac{\theta }{2}\right)|a_j{⟩}_n+\sin \left(\frac{\theta }{2}\right)e^{i\varphi }|b_j{⟩}_n\right).$$

$$|\theta ,\varphi ⟩{⟩}_j=\frac{1}{\sqrt{N!}}(\cos \left(\frac{\theta }{2}\right)a_j^{†}+\sin \left(\frac{\theta }{2}\right)e^{i\varphi }{b}_j^{†}{)}^{N_j}|0⟩.$$

$$H=g{S}^z(J_1^z+J_2^z),$$

Full text

††preprint: APS/123-QED

I Introduction

Quantum phenomena are typically associated with the microscopic world where delicate quantum states are thought to be impossible to generate for all but the smallest number of particles. In the last few decades, macroscopic systems have been created in which large number of particles behave quantum mechanically. Some remarkable examples which have been achieved experimentally are cantilevers being in superposition of oscillation modes O’Connell et al. (2010), non-classical state generation and teleportation of macroscopic atomic ensembles Julsgaard et al. (2001); Riedel et al. (2010); Krauter et al. (2013); Pezzè et al. (2016). On the other hand it is well-known that macroscopic quantum systems often suffer from decoherence exponentially with particle number Schlosshauer (2005); Zurek (2003) and the quantum mechanical properties of the state can only be sustained for very short times. This apparent inconsistency is resolved by understanding that decoherence is a state-dependent process. Some states such as Schrodinger cat states are highly susceptible to decoherence, while coherent states are relatively robust. By using quantum states which are less affected by decoherence, this gives the possibility of realistically realizing macroscopic states that can be used for various applications such as quantum metrology and quantum information Byrnes (2013); Byrnes et al. (2015, 2012); Filip (2013); Ilo-Okeke and Byrnes (2014); Pyrkov and Byrnes (2014a); Tichy et al. (2016); Pyrkov and Byrnes (2014b).

Currently, the types of states used in large ensemble states fall under two basic categories, continuous variable (CV) Furusawa et al. (1998); Braunstein and van Loock (2005) or discrete quantum states. CV has been extensively used to model entanglement and teleportation using degrees of freedom analogous to quadratures of light, and has been extremely successful in demonstrating entanglement generation and teleportation between multiple ensembles Julsgaard et al. (2001); Krauter et al. (2013). Here the transversal spin operators, normally ($J^{z}$ and $J^{y}$) of a $J^{x}$ polarized spin ensemble are treated as canonical position and momentum quadratures of the quantum harmonic oscillator. Deviations from the polarized state on the Bloch sphere, are mapped onto position and momentum variables according to the Holstein-Primakoff (HP) approximation $x=J^{y}/\sqrt{2N}$ and $p=J^{z}/\sqrt{2N}$, where $N$ is the particle number. The approximation as a single quantum harmonic oscillator mode is accurate as long as the quantum states is close to the original $J^{x}$ polarization (Fig. 1(b)). However, for states that deviate from near the $J^{x}$ polarization the approximation breaks down and the full spin degrees of freedom must be considered (Fig. 1(c)). We call such a regime where the HP approximation can no longer be used the ”beyond-CV” regime, as a correct description must involve keeping track of spin, rather than mode operators.

In experiments such as Ref. Julsgaard et al. (2001), the type of entangled state is a two-mode squeezed state under the HP approximation. In this experiment, described schematically in Fig. 1(a), the entanglement is generated by letting two polarization modes of a laser interact with two atomic ensembles via the ac Stark shift Kuzmich and Polzik (2000); Duan et al. (2001, 2000); Cerf et al. (2007). A pair of measurements is then made on the light after interfering the modes. While it is well-established that the procedure creates two-mode squeezed states in the ensembles, given that the ensembles are not genuine quantum harmonic oscillator modes, a question remains of what happens in the light-matter interaction is strong enough such that the HP approximation breaks down. This question is particularly relevant as recently several experiments have reported the generation of non-Gaussian spin ensemble states Gerving et al. (2012); Strobel et al. (2014); Barontini et al. (2015); McConnell et al. (2015); Bohnet et al. (2016); Pezzè et al. (2016). For our purposes we define non-Gaussian states as those with quasiprobability distributions (i.e. the Wigner or Husimi $Q$-distribution) that cannot be described by a Gaussian form. An important question in this context is: what are the characteristics of the states that are generated by the scheme of Fig. 1(a) in the beyond-CV regime?

In this paper, we investigate the quantum states created by a similar protocol to that introduced in Ref. Julsgaard et al. (2001) but without using the HP approximation. The procedure is described in Fig. 1(a), where initially atom-light-atom entanglement is created due to an ac Stark shift interaction. The HP approximation restricts the region of validity to only short interaction times between the atoms and the light. Using our exact time evolution techniques, this allows us to consider long interaction times which can produce states that requires one to use the full spin formalism. The light is then measured, which produces atom-atom entanglement. We analyze the resultant state by examiing the probability distributions, Husimi $Q$-distributions, and entanglement.

This paper is structured as follows. In Sec. II we first review the standard approach using the HP approximation, and derive a procedure that maps between quantum states on the atomic ensembles and the approximated mode representation. In Sec. III we construct a theory of measurement induced entanglement for the scheme shown in Fig. 1(a). Our theory of entanglement generation between two atomic ensembles works in the regime where the HP approximation no longer is valid. In Sec. IV, we obtain approximate expressions by analytic methods to find the relation between the two-mode squeezed state and the qubit entangled state. In Sec. V we plot the probability density and Husimi $Q$-distributions Husimi (1940) of the entangled state. In doing so, we reveal the properties of the entangled state in the beyond-CV regime, where it differs significantly from a two-mode squeezed state. In Sec. VI we analyze the atomic entangled state in terms of the entanglement generation between the ensembles. We finally give our conclusions in Sec. VII.

II Continuous variable approach to light mediated atomic ensemble entanglement

We now review the standard approach of describing atomic ensemble entanglement based on the technique as described in Refs. Kuzmich and Polzik (2000); Duan et al. (2000); Julsgaard et al. (2001); Gross (2012); Hammerer et al. (2010). The approach is based on making a HP approximation to the spin variables of the ensembles, and results in two-mode squeezing interaction. After describing the general framework and results of this approach, we describe what this approximation means in terms of the atomic states such that the approach can be compared to our methods introduced in later sections.

II.1 Standard continuous variables approach

In this section we review the standard theoretical approach of atomic ensemble entanglement based on the HP approximation. In this approach, the discrete collective spin operators are mapped to the continuous canonical position and momentum operators of the harmonic oscillator with unit mass,

[TABLE]

which obey the usual commutation relation $[x,p]=i$. This approximation maps a part of the Bloch sphere around the polarized spin vector in which deviations are discrete into the flat phase space of continuous variables, seen in Fig. 1(b). Intuitively it is clear that this approximation only is valid when the angular deviation is small, as for larger deviations the mapping become increasingly unsatisfactory, as seen (exaggerated) in Fig. 1(c). In this framework, an entangled state is modeled as a two-mode squeezed state in which the variance of the total position and relative momentum tend to zero with increased squeezing. The atom-light interaction in the scheme described in Fig. 1 creates a time evolution in which the states of light and matter are entangled, and a measurement of one reveals information about the other.

In this scheme, an interaction Hamiltonian of the form

[TABLE]

is used, where $S^{z}$ is the Stokes operator describing a light source and $J^{z}$ are the atomic Schwinger boson operators describing the atomic systems. Using the same mapping as in Eq. (1) for both the atomic and photonic operators, we define the canonical continuous variable operators as $x_{\text{at}}$, $p_{\text{at}}$ and $x_{\text{ph}}$, $p_{\text{ph}}$ respectively. As the laser pulse passes through the atomic samples, the photonic operators evolve according to Heisenberg’s equation of motion,

[TABLE]

Assuming the operators remain stationary in time for the duration of the experiment, we get the output modes by simple integration as,

[TABLE]

Similarly, the atomic operators evolve as Kuzmich and Polzik (2000); Duan et al. (2000)

[TABLE]

After the evolution a measurement of $x_{\text{ph}}$ is performed, it is measured and collapses to a constant. According to (3) and the fact that

[TABLE]

from (7) and (8) this means that the quantity $p_{\text{at},1}(t)+p_{\text{at},2}(t)$ also collapses to a constant, provided $x_{\text{ph}}(0)$ is small and the dimensionless interaction parameter $gt/\hbar$ is large. This means that the variance of $p_{\text{at},1}(t)+p_{\text{at},2}(t)$ is small, which results in correlations between these observables. Meanwhile, from (5) and (6) we have

[TABLE]

which means that the correlations between $p_{\text{at},1}(t)$ and $p_{\text{at},2}(t)$ can be produced without changing the correlations between $x_{\text{at},1}(t)$ and $x_{\text{at},2}(t)$.

The procedure is then repeated again but in a different basis in order to generate the correlations between $x_{\text{at},1}$ and $x_{\text{at},2}$. The two atomic systems are rotated around the the $J^{x}$-axis, which transforms the canonical operators for the first system as $x_{\text{at},1}\to-p_{\text{at},1}$ and $p_{\text{at},1}\to x_{\text{at},1}$, while the second system is rotated in the opposite direction, $x_{\text{at},2}\to p_{\text{at},2}$, and $p_{\text{at},2}\to-x_{\text{at},2}$. A second measurement of the light field, $x_{\text{ph}}$ under the same conditions as before fixes the quantity $x_{\text{at},2}-x_{\text{at},1}$ to a constant, and gives a small variance. Similiarly to Eqs. (5)-(8) the evolved quantities $p_{\text{at},1}+p_{\text{at},2}$ and $x_{\text{at},2}-x_{\text{at},1}$ remain conserved through the procedure.

II.2 Mapping between atomic spin states and the HP approximated mode representation

In the continuous variable approximation, one of the total spin operators $J^{x}$ is set to a constant and the other two operators $J^{y},J^{z}$ are treated as approximate position and momentum variables. Typically the effect of the atom and light interaction is derived in the Heisenberg picture where it is relatively simple to derive the entanglement. In this section we describe what this approximation means in the Schrodinger picture, which will be the framework that will be used in our analysis in following sections. We deduce the mapping between the original spin and HP approximated bosonic mode Fock states.

First let us write the total spin operators in terms of bosonic modes, which are conventionally defined as

[TABLE]

These obey commutation relations $[J^{l},J^{m}]=2i\epsilon_{lmn}J^{n}$ where $\epsilon_{lmn}$ is the Levi-Civita antisymmetric tensor. Here, $a^{\dagger}$ and $b^{\dagger}$ denote the creation operators for two hyperfine levels of the atoms in each of the ensembles labeled by $j$. These operators obey bosonic commutation relations $[a,a^{\dagger}]=[b,b^{\dagger}]=1$. To explicitly observe the HP approximation in terms of bosonic modes, we make a SU(2) transformation

[TABLE]

which transforms the spin operators as

[TABLE]

The HP approximation requires that we are in the vicinity of a fully $J^{x}$ polarized state, which is

[TABLE]

Due to the macroscopic population of the $b_{x}$ state, we may take this mode to be a constant $b_{x}\sim\sqrt{N}$. We thus have

[TABLE]

where the position and momentum operators are $x=(a_{x}+a_{x}^{\dagger})/\sqrt{2}$ and $p=(a_{x}-a_{x}^{\dagger})/\sqrt{2}i$ Braunstein and van Loock (2005). We thus see that the HP approximation amounts to setting macroscopically occupied mode to a constant, and treating the remaining mode quantum mechanically.

Given a particular quantum state, we may then transform between the original spin representation and the HP approximated state in the following way. Consider an arbitrary state of the spin ensemble

[TABLE]

where

[TABLE]

are the Fock states for the spins. According to the above procedure, we must set the $b_{x}$ mode to a constant, hence we must first make a change of basis to the $J^{x}$ Fock states defined as

[TABLE]

The state is then written

[TABLE]

where $\psi^{\prime}_{k_{x}}=\sum_{k}\langle k_{x}|k\rangle\psi_{k}$. The matrix elements $\langle k_{x}|k\rangle$ are calculated according to Refs. Tinkham (2003); Thompson (2008). Setting $b_{x}=\sqrt{N}$ in (18) results in an unnormalized state, hence we make the association

[TABLE]

We may thus say that the spin Fock states in the $J^{x}$ basis are the bosonic mode Fock states in the HP approximation.

III Theory of light mediated atomic ensemble entanglement

In this section we introduce our approach of deriving the ensemble-ensemble entanglement in the Schrodinger picture. As shown in Fig. 1(a), linearly polarized coherent light illuminates two ensembles, each in a coherent spin state. A beam splitter then interferes the light, after which the photons are measured. As described in Sec. II.1, the entanglement procedure makes use of two successive measurements of the light field, with the atomic ensembles rotated in between in order to infer correlation in both the $J^{z}$ and $J^{y}$ directions. As the two procedures are identical up to a basis rotation, we restrict our analysis to one of the measurements, which will prove to have rich dynamics as the system is evolved beyond the HP approximation.

The initial state of light is polarized in the $S^{x}$-direction, which is a superposition of left and right circularly polarized light. The initial quantum state of the light can thus be written

[TABLE]

where $\xi$ is the amplitude of the light, and $c^{\dagger},d^{\dagger}$ are the creation operators for left and right hand circularly polarized light. The Stokes operators are defined as

[TABLE]

The initial state of light has the expectation value $\langle S^{x}\rangle=|\xi|^{2}$ but is zero for $\langle S^{y,z}\rangle=0$.

For the atomic ensembles, the initial state is

[TABLE]

where we have defined the coherent spin states as

[TABLE]

Here $|a\rangle_{n},|b\rangle_{n}$ denote the two hyperfine states of the $n$th atom in the ensemble, and $\theta$ and $\phi$ are two arbitrary angles on the Bloch sphere, $0<\theta<\pi$, $0<\phi<2\pi$. The index $j=1,2$ labels the two ensembles, and $N_{j}$ is the number of atoms in each ensemble. In our calculations, it will be convenient to work with the bosonic formulation of coherent spin states, which is valid as long as the wavefunction in the ensemble is symmetric under particle interchange. This will be always true in our case, as the Hamiltonian used to evolve the system and the initial states are symmetric. The spin coherent state in the bosonic formulation is written

[TABLE]

where the bosonic operators $a_{j},b_{j}$ are defined as in Sec. II.2, with the additional $j=1,2$ labels for each ensemble. We note that in Ref. Julsgaard et al. (2001) in fact the ground state of the underlying atoms is a $F=4$ state, and our case above would strictly speaking correspond to $F=1/2$. While these may appear different, another way to view each ensemble is that it is a macroscopic spin with total spin $J=N$, where all the underlying spins are symmetric under particle interchange. In this picture the constituent particles making up the macroscopic spin become irrelevant to the dynamics as long we make observations in the total spin variable. In this sense our results should be also valid for any $F$ making up the ensemble.

With the use of the spin operators we model the interaction between the two atomic ensembles as a non-linear quantum non-destructive (QND) Hamiltonian of the form Kuzmich et al. (2007); Takahashi et al. (1999); Takeuchi et al. (2005)

[TABLE]

where $g$ is the interaction parameter resulting from the ac Stark shift coupling. In terms of physical parameters, this is given by $g=\frac{\sigma\gamma\alpha_{v}}{A(I+\frac{1}{2})\Delta}$ where $\sigma$ is the resonant absorption cross section for an unpolarized photon on an unpolarized atom, $A$ is the cross section area of the light beam, $\gamma$ the spontaneous emission rate from the upper atomic level, $\alpha_{v}$ the vector polarizability and $I$ the value of the nuclear spin Kuzmich and Polzik (2000). The interaction is considered to be a QND if $[H,S^{z}]=0$, which is satisfied in our case.

We evolve the system in time by applying the unitary time evolution operator $U=\exp(-itH/\hbar)$ which yields, with the dimensionless entanglement time $\tau=gt/\hbar$,

[TABLE]

where we have expanded the atomic states in the $J^{z}_{i}$ eigenbasis. The light component of this can be evaluated by noting that a number operator on a state causes a coherent state to pick up a phase $e^{i\theta c^{\dagger}c}e^{\xi c^{\dagger}}|0\rangle=\exp[\xi e^{i\theta}c^{\dagger}]|0\rangle$. This shifts the phase of the optical state with a value dependent on the $J^{z}_{j}$ eigenstate. Now the light and atoms become entangled, giving

[TABLE]

The phases picked up by the interaction are now interfered using a using a 50:50 beam splitter after the light pulse has passed through the two atomic ensembles. The beam splitter transforms the photonic operators as

[TABLE]

which yields

[TABLE]

where $x=\tau(2k_{1}+2k_{2}-(N_{1}+N_{2}))$. The last step is to project the above state on the photonic number states $|n_{c},n_{d}\rangle$, with measurement outcomes $n_{c}$ and $n_{d}$ respectively. We note here that depending on whether a photon resolving measurement is made, $n_{c}$ and $n_{d}$ may or may not be explicitly known. Regardless of whether this is known, we shall see that entanglement will be produced between the atomic ensembles, on a shot-to-shot level.

The state after projecting on the photonic number states is

[TABLE]

where

[TABLE]

and

[TABLE]

is a normalization factor as a measurement of the light component was made, and irrelevant global phase factors were dropped. Eq. (31) is the core result of this paper, the latter sections are devoted to examining its properties. For zero interaction time $\tau=0$, $A_{n_{c},n_{d}}(k_{1},k_{2})$ is constant in $k_{1}$ and $k_{2}$ and the state of the atoms are unaffected by the measurement. Evolving in $\tau$ gives entanglement between the two ensembles due to the correlations in $A_{n_{c},n_{d}}(k_{1},k_{2})$. In the next section we analyze this function to see what kind of correlations are present due to the atom-light interaction.

IV Analytic approximation of probability densities

As Eq. (31) is the wavefunction after the interaction with the light, in principle it is possible to derive all physical quantities based on this. However, due to the complicated nature of the coefficient $A_{n_{c}n_{d}}(k_{1},k_{2})$ it is not entirely obvious what kind of entangled state this is and hence it is beneficial to make some analytical approximations. For small light-ensemble interaction times, we expect the atomic entangled state to behave like a two-mode squeezed state, as the HP approach should be a good approximation in this regime. First, the binomial coefficient can with good approximation be treated as Gaussian when $N\gg 1$,

[TABLE]

Secondly, using Stirling approximation we may write for $N_{p}\gg 1$ the combination of the binomial and trigonometric functions as a Gaussian function Ilo-Okeke and Byrnes (2016),

[TABLE]

with $x=2\tau(k_{1}+k_{2}-N)$. Here, the measured photon numbers in each mode are $n_{c}$ and $n_{d}$ and thus the total photon number is $n_{c}+n_{d}=N_{p}$. The maximum of the Gaussian occurs at

[TABLE]

which has in general multiple solutions for $k_{1},k_{2}$, due to the oscillatory nature of the sine function. As the beam splitter in Fig. 1(a) is a 50:50 beam splitter, on average we expect that the number of photons in each detector is equal $n_{c}\approx n_{d}\approx N_{p}/2$. Restricting our analysis to small times $\tau\sim 1/N$, we may linearize the trigonometric functions around $x\ll 1$, and obtain the simplified expression for the probability distribution,

[TABLE]

The first factor in this expression states that atom Fock states occur with averages $\langle k_{1,2}\rangle=N/2$ and with standard deviation $\sqrt{N}/2$. This is the same distribution as the original spin coherent states (23) that the atoms are initially prepared in. The second factor gives a correlation between the two atomic modes which results in entanglement. The strength of the correlation is dependent on the square of the dimensionless entanglement time and the number of photons in the interaction.

To confirm our analytical prediction, we numerically plot the probability distribution of the wavefunction in Eq. (31) for different entanglement times in Fig. 2. When $\tau=0$, there is no correlation between the two modes and the probability distribution is a Gaussian distribution centered around $k_{1}=k_{2}=N/2$. Evolving the interaction to times of order $\tau\sim 1/N$, the state develops correlations between Fock states of the ensembles, in agreement with (37). This state should be consistent with a two-mode squeezed state, according to the arguments of Sec. II.1. To verify this, let us write the wavefunction of the two-mode squeezed state

[TABLE]

where the $|n\rangle$ are bosonic mode Fock states. To compare the probability distributions of the states, we employ the mapping procedure as discussed in Sec. II.2, where the mode Fock states are mapped back onto the spin Fock states in the $J^{x}$ basis. We thus expect that the two-mode squeezed state has a state in terms of the spin states according to

[TABLE]

where we have restricted the sum to $N$ as now these are now Fock states for the spins, and ${\cal N_{\text{TMS}}}$ is a suitable normalization constant. In the limit $N\gg 1$, it is possible to evaluate the sum over $k$ analytically (see Appendix), and we obtain

[TABLE]

We see an immediate similarity of the mapped two-mode squeezed state probability distribution to (37). As with the case of the atomic entangled state, there is a Gaussian envelope from the first exponential factor with $\langle k_{1,2}\rangle=N/2$ and standard deviation $\sqrt{N}/2$. The second exponential induces a correlation between the two modes $k_{1}$ and $k_{2}$ with squeezing parameter $r$. Increasing the squeezing parameter $r$ causes the coefficient of the second Gaussian $\sinh(2r)$ to increase exponentially, narrowing the distribution in the direction $k_{1}-k_{2}-N=0$ while extending it in the opposite direction. Atomic ensembles have a fixed atom number which restricts the range of $k_{1,2}$, while bosonic modes do not have this restriction. This means that the mode Fock state occupation probability extends to infinity for the HP approximated case. Inspecting the two probability density functions, we can infer a relation between the entanglement time $\tau$ and the squeezing parameter,

[TABLE]

For small squeezing $r\ll 1$, we may approximate $\sinh(2r)\approx 2r$, hence we obtain $2N_{p}N\tau^{2}\leftrightarrow r$ in this regime. In Fig. 2 and Fig. 3 we compare the probabilities of the atomic states (31) after the interaction with the light with the mapped two-mode squeezed state. We see in general good agreement between the two distributions, with the expected correlations developing with squeezing parameter $r$.

Further increasing the entanglement time $\tau$, the Fock state correlations become stronger, which the distribution narrowing as would be expected of (37). For interaction times of order $\tau\sim 1$, multiple probability distribution peaks appear indicating non-Gaussian behavior (Fig. 2(c)(d)), due to the periodicity of the trigonometric functions in (36). The boundary between the linear regime where only one and many correlation peaks occurs is the time $\tau\sim 1/\sqrt{N}$, where weak secondary peaks start to appear in Fig. 2(c). The corresponding exact probability distribution is shown in Fig. 4(a), showing good agreement to Fig. 2(c) After these times the two probability distributions diverge, and the atomic entangled state entangles $k$-states other than $k_{1}+k_{2}=N$, which is apparent from Fig. 2(c)(d). These correlations beyond the primary peak emerge when there are more than one solution for the maximum of (35). The multiple peaks signal the breakdown of the HP approximation, where the correlations cannot be treated linearly.

One difference between the approximated distribution and the atomic entangled state is that probability distributions of the former tend to broaden as the number correlations become stronger, while the latter distributions stay within a fixed envelope. This occurs because in the infinitely squeezed limit the mapped two-mode squeezed state (39) is

[TABLE]

The probability distribution is thus a perfectly Fock number correlated state with a uniform distribution from $k=0$ to $N$. Such a state never truly occurs in the atomic ensemble wavefunction, due to the development of multiple correlation peaks within a finite envelope probability. The closest distribution to this would be that shown in Fig. 2(b), where there is one correlation peak and the times are still within the HP approximated regime. From (41) we can estimate the maximum squeezing that is possible

[TABLE]

This suggests that it is advantageous to use very bright light in comparison to the ensemble particle number in order to achieve a high squeezing. We note that (43) should be viewed as a theoretical upper bound to the squeezing and in practice there will be other factors such as decoherence which will reduce the maximum possible squeezing further.

At even longer times, the lines become increasingly dense until $\tau=\frac{\pi}{8}$, where the probability becomes a checkerboard pattern, where every second $k$-state has a non-zero probability, as seen in Fig. 4. Inserting this time, $\tau=\frac{\pi}{8}$ into (31), while assuming $n_{c}=n_{d}$, a sum of trigonometric functions, with $m=k_{1}+k_{2}-N+1$ is obtained,

[TABLE]

Depending on the parity of $n_{c}$, this evaluates in two different ways. For even $n_{c}$, \cos^{n_{c}}\big{(}\frac{m\pi}{2}\big{)}=\frac{1}{2}(1+(-1)^{m}), while for odd $n_{c}$, \cos^{n_{c}}\big{(}\frac{m\pi}{2}\big{)}=\frac{1}{2}(i^{m}+(-i)^{m}). We then have, for odd $n_{c}$,

[TABLE]

Dropping overall constant factors and using Eq. (24), we can identify the sums as coherent spin states at distinct positions on the Bloch sphere. The odd state then becomes

[TABLE]

while the case for even $n_{c}$ results in a similar state,

[TABLE]

Now we see that for this particular time, assuming the measured photon numbers in each mode are equal, we obtain Bell states consisting of states at opposite sides of the Bloch sphere. Since the superposition of states involves spin coherent states at opposite sides of the Bloch sphere, we call the states (46) and (47) ”Schrodinger cat-Bell states”. However, the axis in which the states are polarized depends on the parity of the $n_{c}$. For odd values, the state end up in $\pm J^{y}$, while for even numbers, the direction is switched to $\pm J^{x}$. A similar form of entangled states was obtained for a pure $J^{z}_{1}J^{z}_{2}$ interaction in Ref. Byrnes (2013). At a time $\tau=\frac{\pi}{4}$ the atomic states become disentangled again, and the states are unchanged from the initial state $|\psi_{\text{atoms}}(t=0)\rangle$.

V Husimi $Q$-distributions

In the previous section we saw that the light mediated entanglement produced a state with a probability distribution consistent with a two-mode squeezed state for small interaction times, and a more complex correlation structure for longer times. In order to understand the nature of the correlations we analyze the atomic entangled state in terms of coherent spin state Husimi $Q$-distributions, defined as

[TABLE]

As the state that we wish to analyze has two ensembles, the $Q$-function has four degrees of freedom $\theta_{1},\phi_{1},\theta_{2},\phi_{2}$, and cannot be visualized directly. For this reason we consider the marginal and joint $Q$-functions, where one of the ensembles is traced over in the former case, and the conditional distribution is examined in the latter. The resulting $Q$-distributions then become a function of two variables which can be visualized.

V.1 Marginal $Q$-functions

The marginal $Q$-function is obtained by tracing out one of the ensembles of the atomic entangled state

[TABLE]

and finding the overlap with coherent spin states as defined in (25). We obtain

[TABLE]

In Fig. 5 we see the behavior of the marginal $Q$-function for different interaction times. The $Q$-function starts as a symmetric Gaussian at $\tau=0$, as the initial state is simply a coherent spin state and is independent of the second ensemble Byrnes et al. (2012)

[TABLE]

For short times $\tau=1/N$, the $Q$-function starts to broaden in the $\phi$ direction, due to the $J^{z}_{1}+J^{z}_{2}$ Hamiltonian which rotates the spins around the $z$ axis of the Bloch sphere. The partial trace has the effect of collapsing all the various terms in the entangled state, creating a broadening effect, in a similar way to Ref. Byrnes (2013). Up until the characteristic time $\tau\lesssim 1/\sqrt{N}$, the coherent states are all distributed around the original point in the Bloch sphere to a high degree, due to the fact that $A(k_{1},k_{2})$ in (32) varies little for different $k$ in the sum. Further evolving in $\tau$ creates Gaussians centered at different positions around the equator of the Bloch sphere, in analogy to that seen in Ref. Byrnes (2013). At $\tau=\pi/8$ there are two Gaussian distributions at $\phi=0,\pi$ and $\theta=\pi/2$, which we attribute to the Schrodinger cat entangled state (47). Finally the $Q$-functions are periodic in time $\tau=\pi/4$, as would be expected from (25).

Let us now compare the $Q$-functions to a genuine two-mode squeezed state to see the degree of agreement. The marginal $Q$-function (in the space of bosonic modes) of the two-mode squeezed state is

[TABLE]

where the coherent state is

[TABLE]

and

[TABLE]

This gives the $Q$-function

[TABLE]

which is a Gaussian function of radius $\sim 1/(1-\tanh^{2}(r))$. We compare this by applying the spin-mode mapping from section II.2 to calculate the $Q$-function (for bosonic modes) of the mapped atomic entangled state. The first step of the mapping is to transform the state into the $J^{x}$ basis, giving

[TABLE]

We then interpret the $|k_{x}\rangle$ states to be the bosonic mode Fock state, which allows us to evaluate the $Q$-function to be

[TABLE]

This distribution is defined in the phase space of bosonic modes as a function of the complex parameter $\alpha$.

The comparison between the two distributions is shown in Fig. 6. We see that for the genuine two-mode squeezed state, the marginal $Q$-distribution always remains a symmetric Gaussian, with its radius increasing with squeezing parameter $r$. The reason for this behavior can be seen from the density matrix (54) which shows that the state is a completely mixed state in the mode Fock states. For $r=0$, the state is a vacuum state, and smoothly evolves towards a mixed state involving all mode Fock state numbers. In contrast, the dominant behavior of the atomic entangled state is to broaden along the equator, but not in the longitudinal direction. This is as expected as in our expression (31) we only perform the first entangling step producing correlations in the $\phi$ direction but not $\theta$. A second measurement would produce the correlations in the $\theta$ direction, with a corresponding broadening in this direction. The spin and mode $Q$-functions show a general similarity (Figs. 5(a)(b) and 6(a)(b)), particularly in the short time region where the HP approximation is valid.

V.2 Ensemble-ensemble correlations

One of the most striking features of the two-mode squeezed state is the presence of non-local correlations between the ensembles, make it useful for tasks such as quantum teleportation. For an infinitely squeezed two-mode squeezed state (42), projecting on a $|k\rangle$ or $|k_{x}\rangle$ state on one mode collapses the other mode into the same state. In this section we examine the correlations between the modes due to the entangling procedure, and compare them to a genuine two-mode squeezing operation.

First let us analyze what is expected for the $Q$-function of a genuine two-mode squeezed state. Evaluating the $Q$-function (for bosonic modes) with respect to (38) gives

[TABLE]

In the limit of strong squeezing $r\rightarrow\infty$, we have $\tanh(r)\rightarrow 1$, and only the first exponential factor is present. If we now set $\alpha_{2}$ to a constant and consider the $Q$-distribution on mode 1, we observe that it is centered around $\alpha_{2}^{*}$. In terms of position and momentum, this corresponds to

[TABLE]

i.e., the positions are correlated and momenta are anti-correlated. For finite squeezing, the additional Gaussian factor in (58) suppresses the correlation, and the relations (59) tend to saturate at finite values (see Fig. 7(a)(b)).

With this in mind, let us examine the $Q$-functions of the atomic entangled states. We first examine the spin $Q$-functions, calculating the overlap of the state with respect to coherent spin states (48), for which we obtain

[TABLE]

The a entangling times are shown in Fig. 7(c)(d)(e)(f), for various choices of interaction times $\tau$ and coherent spin state positions $\theta_{2},\phi_{2}$. We choose the coherent state centers on ensemble 2 according to where the marginal distributions have a significant probability. Choosing $\theta_{2},\phi_{2}$ that are outside the marginal distribution corresponds to very low probability events, giving distributions that are more difficult to interpret. In Fig. 7(c), we choose a coherent state center displaced along the equator in the positive direction. We see that as expected the coherent spin state follows the chosen angle $\phi_{2}$ by the same amount. This is the same behavior as would be expected of a two-mode squeezed state. In Fig. 7(d), we vary the polar angle $\theta_{2}$ instead. For small times $\tau=1/N$, we see that the distributions shows a weak anit-correlation in the angle. The anti-correlation is weaker than in the $\phi$ direction, and does not preserve the symmetrical nature of the $Q$-function. This is again due to the fact that we only include the first measurement in our analysis, which produces the correlations in the $\phi$ direction. It is however interesting that a despite not putting in the correlation at all, some anti-correlation is present. For longer times $\tau=1/\sqrt{N}$, the $\theta$ correlation becomes weaker, although the $\phi$ correlation remains in place.

To better quantify the correlation and anti-correlation between the two ensembles, we take the spin expectation values on ensemble 1 after projecting a coherent spin state on ensemble 2. Denoting the projector on the second ensemble as $P(\theta_{2},\phi_{2})=|\theta_{2},\phi_{2}\rangle\rangle\langle\langle\theta_{2},\phi_{2}|$, the expectation values are calculated as

[TABLE]

where we require a normalization of the state as a projection is being made. We obtain expressions for the average spins as

[TABLE]

and

[TABLE]

where ${\cal N}_{\theta\phi}$ is the normalization factor, common for both expectation values,

[TABLE]

The expectation values are shown in Fig. 8 for two different entanglement times, $\tau=\frac{1}{N}$ and $\tau=\frac{1}{2\sqrt{N}}$. In Fig. 8(b) we show the correlation with $\phi_{2}$. As expected we see a positive correlation between the coherent spin state $\phi_{2}$ and the average spin $\langle J^{y}_{1}\rangle$. According to (59), for a genuine two-mode squeezed state we should observe under the HP approximation, Eq. (1),

[TABLE]

where the expectation value involves the projection $P(\theta_{2},\phi_{2})$. For both the interaction times, the $\langle J^{y}_{1}\rangle$ expectations follow the expected sine curve. The $\langle J^{z}_{1}\rangle=0$ indicating that the state remains on the equator of the Bloch sphere at all times. The longer times start to deviate from the genuine two-mode squeezed state, due to the non-Gaussian characteristics that are developing in the state. For the $\theta$ correlations (see Fig. 8(a)), we see the anti-correlations are generally weaker as expected. While in the genuine CV regime one would not expect to see correlations at all (as only the first measurement is included), we nevertheless see a positive slope in which an increase in $\theta$ from the ground state $\theta=\frac{\pi}{2}$, corresponding to a negative change in the mapped momentum operator, $p$ results in a positive expectation value. The anti-correlations diminish with larger interaction times, hence is generally a weaker effect. Similarly to the $\phi$ correlations, the $\langle J^{y}\rangle$ is equal to zero for all entanglement times, indicating there is no correlation with $\phi$, as expected.

VI Entanglement

In the previous section, we showed that as the interaction times grow longer, the states evolve into complex states with non-Gaussian properties. For small interaction times the state has the characteristics of a two-mode squeezed state, but for longer times the state deviates from this and the expected correlations degrade. One may be tempted to conclude from this observation that entanglement is degrading in the system. This is in fact untrue as the two-mode squeezing interaction is only an approximation at short times, and at longer times the entangling scheme produces a different type of entanglement at longer times. Due to the large number of degrees of freedom in an ensemble, it is possible to have many different types of entangled state, which are not necessarily simply described as a two-mode squeezed state. This is a similar situation to Ref. Byrnes (2013), where the system evolves according to a $J^{z}_{1}J^{z}_{2}$ interaction, and a variety of different types of entangled states are produced.

To show this explicitly, we plot the entanglement as quantified by the von Neumann entropy

[TABLE]

where the density matrix is given in (49). The results are shown in Fig. 9. We see that the entanglement shows a remarkably complex behavior with large dips at certain times. In Fig. 9(a) the measured photon numbers are equal, and the entanglement is symmetric around the characteristic time $\tau=\pi/8$. For $N=1$ (i.e. a single atom in each ensemble) we see that we recover maximal entanglement at $\tau=\pi/8$, where the state has the form of a Bell state. The saturation of the entanglement at around half the maximal entanglement is exactly the same behavior as that seen in Ref. Byrnes (2013), which arose due to the binomial factors present in the wavefunction, which is also present in our case (31). The binomial factors create an uneven superposition of states which only involve spin Fock states that are centered around $k_{1},k_{2}=N/2$ with width $\sim O(\sqrt{N})$. We attribute the saturation here to the same origin. We also plot the entanglement of an ideal two-mode squeezed state (39), using the mapping (41) to convert the squeezing parameter $r$ to interaction time. For the ideal two-mode squeezed state, the entanglement increases towards a maximally entangled state, as can be directly seen from (39) where there is an equal superposition of pairs of spin Fock states. Such a state is never attained in our protocol due to the binomial factors as discussed above.

As the measured photon numbers are changed, we see in Fig. 9(b) that the entanglement follows overall the similar behavior, but different amplitudes of the dips in the entanglement. As the change in photon numbers does not change the arguments in Eq. (32), the dips appear at the same times. This behavior is remarkably similar of the fractal ”devil’s crevasse” entanglement produced by a $J^{z}_{1}J^{z}_{2}$ interaction Byrnes (2013). Although in this case the interaction is not precisely the same, the Hamiltonian (26) shares the same form and we attribute the similar entanglement to similar resonances of coherent state as seen in Ref. Byrnes (2013). The complex structure of dips in the entanglement occur only after the characteristic time $\tau=1/2\sqrt{N}$. Prior to this time there is a monotonic increase of the entanglement, with little variation on the photon number outcomes. We thus expect that for times $\tau>1/2\sqrt{N}$ there will be a strong dependence of the state on the obtained photon number. Statistical averaging of states at this time will likely produce states with lower purity, which will degrade the entanglement, if no postselection is used. In contrast, for earlier times $\tau<1/2\sqrt{N}$, there should be less dependence on the photon number. In the limit of very short interaction times $\tau\sim 1/N$ the scheme reduces that to the CV approximated case, which should have very little photon number dependence to the entangled state that is created. Such a behavior is consistent with the results found in Ref. Byrnes (2013), where the states where robust in the presence of decoherence for entangling times in the region $\tau\sim 1/N$, but sensitive for longer times $\tau\sim O(1)$. The crossover between decoherence robust and sensitive regimes occurred at times $\tau\sim 1/\sqrt{N}$, which is consistent with the results we observe here.

VII Summary and conclusion

We have presented a model and analyzed the quantum state for the light-mediated entangling procedure as described in Refs. Julsgaard et al. (2001); Kuzmich and Polzik (2000); Duan et al. (2001, 2000); Cerf et al. (2007) with the modification that we restrict us to modeling interaction in one direction. The main benefit of our approach as compared to previous theoretical descriptions is that it does not make the HP approximation, where one of the spin degrees of freedom is treated as a c-number. This allows us to investigate the beyond-CV regime, where the state starts to develop non-Gaussian characteristics. We have investigated this state from the point of view of the Fock number probability distributions, the marginal and joint $Q$-function distributions, and entanglement between the ensembles. From the simulations we find that the entangled state behaves like a two-mode squeezed state for short entanglement times of $\tau\sim 1/N$. Further evolving in time, new behavior is seen where the entangled coherent spin states become non-Gaussian in the marginal $Q$-function and wrap around the equator of the Bloch sphere. As the system enters the non-Gaussian regime, the correlations tend to degrade, where the simple CV correlations between coherent states breaks down. Additionally, we have found an analytical relation between the squeezing parameter when using the HP approximation and the entanglement time of the atomic entangled state, which allows us to compare the genuine two-mode squeezed state to the atomic state. This gives a theoretical maximum squeezing that is possible, defined as where the HP approximation starts to break down.

The diminishing of the two-mode squeezed correlations is not due to additional effects such as decoherence entering in the calculation – our calculations do not take these into account and thus may be considered to be the “ideal” behavior. This is evident by looking at the entanglement between the ensembles as shown in Fig. 9. In fact, the regime where the HP approximation is valid $\tau\sim 1/N$ in fact relatively has small amounts of entanglement compared to longer interaction times. Thus entanglement is certainly present at all times, but it is not a two-mode squeezed state, but a more complex type. Due to the large number of degrees of freedom of an ensemble, there are many different types of entanglement that are possible, in addition to a two-mode squeezed state. These have a non-Gaussian nature, as shown by the $Q$-functions, and potentially have large amounts of entanglement in them. At particular times such as $\tau=\pi/8$ it is possible to write exactly as an entangled Schrodinger cat-Bell states (46) and (47). Such deviations from the HP approximation can also potentially occur in other experimental configurations than those considered here, for example those in Refs. Peise et al. (2016, 2015) where two-mode squeezed states are induced by atomic interactions between hyperfine states. The degree of deviation depends upon the degree of depopulation of the mode which is treated classically.

In this paper we did not include the effects of decoherence, as we found that the pure state properties are already rather rich, and have not been analyzed fully to our knowledge. In a realistic setting, the main sources of decoherence would be photon and atomic losses, which we expect to degrade the entanglement as observed in Ref. Byrnes (2013). If the same behavior is observed as the $J^{z}_{1}J^{z}_{2}$ interaction, then for times $\tau>1/\sqrt{N}$ decoherence effects are greatly enhanced due to the generation of Schrodinger cat-like states, which are known to be highly susceptible to decoherence. This scenario is quite likely due to the similar nature of the interaction Hamiltonian $H\propto S^{z}(J^{z}_{1}+J^{z}_{2})$. This suggests that there is a window of interaction times $1/N\lesssim\tau\lesssim 1/\sqrt{N}$ which is beyond the CV regime and shows non-Gaussian characteristics, yet is still stable under decoherence. One requirement would be that the entanglement state between the atoms would be relatively invariant to the photon counting outcomes $n_{c},n_{d}$. To our knowledge this has yet to be observed experimentally, and would display rich physics with potential applications to quantum metrology and computing Gerving et al. (2012); Strobel et al. (2014); Byrnes et al. (2015, 2012).

VIII acknowledgements

O. P. and T. B. thank Qiongyi He and Yumang Jing for discussions. 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 (Grant No. 61571301), the Thousand Talents Program for Distinguished Young Scholars (Grant No. D1210036A), and the NSFC Research Fund for International Young Scientists (Grant No. 11650110425).

Appendix A Derivation of two-mode squeezed probability distribution

Here we explain how the probability distribution of (39) can be evaluated to (40). In the region where the system is strongly polarized, near the extreme values $k=0,N$, the transformation elements $\langle k_{x}|k\rangle$ can be well approximated as harmonic oscillator functions, as described in Ref. Rowe et al. (2001),

[TABLE]

In this transformation, $m=\frac{N}{2}-|\frac{N}{2}-k^{x}|$, $x_{1,2}=-\frac{1}{\sqrt{2N}}(2k_{1,2}-N)$ and $H_{m}(x)$ is the $m$th Hermite polynomial. Armed with this we change basis from $k_{x}$ to $k$, and obtain,

[TABLE]

where $\lambda=\tanh(r)$. Using the completeness relation for Hermite functions, the sum over $k_{x}$ can be expressed in terms of Mehler kernels Viskov (2008), which states that

[TABLE]

Using the above definitions for $x_{1}$ and $x_{2}$, and denoting $\alpha=\frac{1-\lambda}{1+\lambda}=e^{-2r}$, the two-mode squeezed state probability density in number representation, $p(k_{1},k_{2})=|\langle k_{1}|\langle k_{2}|\text{TMS (mode)}\rangle|^{2}$ can be expressed in terms of Gaussian functions,

[TABLE]

which after rearranging gives

[TABLE]

Some elementary manipulations give (40).