Mesoscopic entanglement through central–potential interactions

The generation and detection of entanglement between mesoscopic systems would have major fundamental and applicative implications. In this work, we demonstrate the utility of continuous variable tools to evaluate the Gaussian entanglement arising between two homogeneous levitated nanobeads interacting through a central potential. We compute the entanglement for the steady state and determine the measurement precision required to detect the entanglement in the laboratory.


Introduction
The mastery over levitated optomechanical systems in the laboratory is becoming increasingly refined. With advances in cooling a small levitated sphere to the ground-state [1][2][3][4][5] and in the preparation of squeezed states [6], compounded by the ability to control and implant charges into levitated nanobeads [7], setups are reaching unprecedented levels of control. Furthermore, optomechanical systems have shown significant potential for sensing applications [8,9], especially with regards to measuring gravitational parameters [10][11][12]. As a result, the engineering of schemes to entangle multiple levitated oscillators and, more broadly, mesoscopic systems, is being identified as a major medium-term milestone of the field [13]. In fact, the entanglement for a number of mesoscopic systems has already been demonstrated experimentally [14][15][16], although a more thorough and exhaustive understanding of * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
the conditions under which entanglement may be generated is required to move on to applications.
A key property of entanglement is that it acts as an unambiguous hallmark of non-classicality. As a result, entanglement between mesoscopic systems can aid the quest of mapping the transition from the quantum to the classical scale [17,18]. Furthermore, recent proposals concerning the fundamental nature of gravity have considered quantum entanglement as generated by a Newtonian potential between two massive quantum systems [19,20]. The Newtonian potential can be seen as an effective potential that arises form an underlying fully quantum field theory. Successfully detecting gravitational entanglement would indeed be a strong indication of the quantisation of gravity [21] and, in general, the detection of entanglement as generated by gravity has significant ramifications, in the attempt to feed theories of quantum gravity with new empirical evidence at low energies [22].
The experimental setup envisioned in [19] is challenging since it requires the generation of highly localised spatial superpositions with large separation for two levitated systems. Gaussian states, on the other hand, are more well-understood, and they can be straight-forwardly prepared in the laboratory. As such, we wish to explore whether the questions posed in [19] can be tested with Gaussian states only. Furthermore, to detect entanglement from the Newtonian potential, it is reasonable to first consider an analogous case with a stronger potential, such as the Coulomb potential.
Our main question in this work is therefore: is it possible to detect entanglement between two mesosopic systems as generated by a central potential interaction with only Gaussian resources? To address this question, we consider a fundamental or effective central potential of the form 1/r n , for integer n and where r = |r 1 − r 2 | with position vectors r 1 and r 2 acting between two spherically symmetric quantum systems. Crucially, we ask whether minimal initial state preparation, such as squeezing, as opposed to severe requirements of preparing highly non-Gaussian states (as, for example, proposed in [19]) is enough to generate detectable entanglement, and whether the witnessing of the generated entanglement is possible simply by measuring position-momentum correlations. We then quantify the leading-order contribution to entanglement within the continuous-variable (CV) framework for both the dynamical generation of entanglement from initial squeezed states of the interacting oscillators and for the system's steady-state in the presence of noise.

Dynamics
We begin by considering two spheres trapped next to each other as per figure 1. The non-interacting system Hamiltonian H 0 that describes the harmonic motion for both spheres in the local trap is given bŷ where m 1 and m 2 are the system masses, ω 1 and ω 2 are the respective mechanical trapping frequencies for each sphere, andx i andp i with i = 1, 2 are the position and momentum operators for system 1 and 2 respectively. We now consider a generic central potential of the form α/|r 1 − r 2 | n , where α is the coupling constant and r 1 and r 2 are position vectors. In all cases considered here, this potential is the lowest-order approximation to a description that includes the full field theory. We proceed to derive the interaction from the generic central potential to second order. The quadratic terms will capture the largest contribution to the central-potential entanglement and retain the quadratic interaction, which maps input Gaussian states to output Gaussian states. The last effect allows us to model the system completely within the covariance matrix formalism [23].
To derive the Hamiltonian interaction term, we assume that the movement of the spheres is constrained in all but the x-direction, which is the axis along which the two systems are trapped. We consider small perturbations to r 1 and r 2 , such that r 1 → r 1 − x 1 and r 2 → r 2 − x 2 , with x 1 r 1 and x 2 r 2 . By denoting r = r 1 − r 2 and Δx = x 1 − x 2 , we Taylor-expand the interaction to second order in Δx to find: where we have left out the dimensional prefactor α and ignored all terms of order O (Δx) 3 . We then quantise the positions of the two masses around their equilibrium positions, which are taken to be a distance r apart from each other, by promoting the position coordinates to operators: x i →x i . For gravity, this step incorporates the assumption that gravity is a quantum force that can path-entangle two quantum systems [19].
To arrive at the entangling interaction term, we discard the constant term, which is the first term on the left-hand side in equation (2), as its only contribution is a static energy shift. Secondly, since a displacement term in a quadratic Hamiltonian does not affect the entanglement [23], we can also discard the term linear in position, which is the second term on the lefthand side in equation (2). The final term, however, contains a mixing of the position operators in the formx 1x2 , which will generate entanglement between the two states. The interaction term in the Hamiltonian thus becomeŝ The full Hamiltonian therefore readŝ SinceĤ is quadratic in the canonical operators, all initial Gaussian states exclusively evolve into Gaussian states. Furthermore, all Gaussian states are uniquely defined by their first and second moments, which allows us to model this system completely within the covariance matrix framework [23]. We introduce the 4 × 4 two-mode covariance matrix σ, which consists of all second moments of the Gaussian stateρ G (t). It is defined as where the vector of operators is given byr = (x 1 ,p 1 ,x 2 ,p 2 ) T , and where the bracket {·, ·} denotes the symmetrised outer product, in the sense that {r,r T } =rr T + (rr T ) T [23]. The evolution of the second moments under the Hamiltonian in equation (4) can be encoded through the symplectic matrix S = e ΩHt/ , where H is the Hamiltonian matrix, defined asĤ = 1 2r T Hr +r T r, ( 6 ) and where Ω is the symplectic form defined in this basis as for a total of n modes. For a bipartite system like the one considered here, n = 2, which means that Ω, σ, H, and S are all 4 × 4 matrices. The covariance matrix σ then evolves as σ(t) = Sσ 0 S T , where σ 0 encodes the second moments of the initial state.
To determine the Hamiltonian matrix H that arises from equation (4), we define the dimensionless operatorsx i andp i aŝ For notational simplicity, we then assume that m 1 = m 2 = m, and ω 1 = ω 2 = ω m . The Hamiltonian in equation (4) therefore becomeŝ In what follows, we will rescale the laboratory-time t by ω m to obtain the dimensionless time parameter τ = ω m t.
We later consider open-system dynamics where κ denotes a mechanical decoherence rate. Here, we also rescale κ toκ = κ/ω m . This yields the Hamiltonian matrixH = H/( ω m ): where for convenience we have defined the dimensionless couplingα and where H 0 = 4 is the 4 × 4 identity matrix that governs the free evolution and H (i) I denotes the Hamiltonian matrices responsible for the interaction. They are given by where specifically H (2) I will generate the entanglement. We note that H (2) I is of the form of two-mode squeezing, which implies that the corresponding closed system will not display periodic behaviour. With these rescaled quantities, the evolution is now encoded as S = e ΩHτ .

Computing the entanglement
To compute the entanglement that arises from the central-potential interaction, we make use of the logarithmic negativity [24][25][26], a well-known monotone that quantifies the degree of violation of the positive-partial-transpose (PPT) criterion [27,28]. The latter is inspired by the by the fact that, given a separable stateρ =ρ A ⊗ρ B , the partial transpose with respect to one of the subsystems leaves the state with positive eigenvalues:ρ Tp 0. Hence, should we find thatρ Tp < 0, the state is entangled. Notice that this criterion turns out to be necessary and sufficient for Gaussian states [29,30].
The PPT criterion in the CV framework can be explicitly computed by dividing σ into submatrices σ A , σ B and σ AB as such [31]: We define the symplectic invariant quantity Δ = det σ A + det σ B + 2 det σ AB . In this basis, the partial transpose is equivalent to settingp i → −p i for one subsystem, which implies The PPT criterion for two-mode Gaussian states can thus be compactly expressed as det σ −Δ + 1 0 [23]. When this equality is violated, the state is entangled. For bipartite states, this is a necessary and sufficient condition for entanglement [29].
To quantify the entanglement, we calculate the logarithmic negativity E N of the state, defined by E N (σ) = max(0, −log 2ν − ), whereν ∓ are the symplectic eigenvalues of the state, defined for bipartite systems as In this work, we consider a two-mode mechanically squeezed state σ S = diag(z, z −1 , z −1 , z) as the initial state. Note that the squeezing occurs in the opposite quadrature, which serves to increase the entanglement from this particular interaction. When z = 1, the state is a coherent state with σ 0 = diag(1, 1, 1, 1). There are many ways to mechanically squeeze the system to produce σ S , including controlling the trap frequency or using a Duffing non-linearity [32]. In addition, thermal optomechanical squeezing has been experimentally realised in [6].
For the system to be thermally stable (i.e. for the Hamiltonian to be bounded from below), we require that the rescaled Hamiltonian matrixH in equation (9) satisfiesH > 0, which means that the eigenvalues ofH must be positive [23]. ForH in equation (9), we find that the eigenvalues λ i for i = 1, 2, 3, 4 are given by λ 1 = λ 2 = λ 3 = 1, and λ 4 = 1 + 2α. This means that for a positive (repulsive) couplingα, there is no limit to the interaction strength that will unbound the Hamiltonian. We are therefore able to obtain considerable amounts of entanglement for a repulsive potential. Whenα is negative (attractive), however, the maximum strength of the interaction is limited tõ α > −1/2. Since all three potentials considered in this work can be attractive (withα < 0), and because strong repulsive potentials pose the risk of eventually losing the trapped beads, we will mainly explore the values that can be obtained for α ∼ −0.4, which is close to the limiting value.
We now proceed to compute the entanglement E N . While an analytic expression for E N is available, it is too long and cumbersome to reproduce here. Instead, we plot the entanglement E N for different values ofα in figure 2. In figure 2(a) we plot E N for coherent states with z = 1 as a function of time τ for various values ofα. The stronger the coupling, the more entanglement is generated. Next, in figure 2(b) we plot the same general dynamics forα = −0.4 but for different z. We

Open system dynamics
All systems are subject to environmental noise that generally degrades the entanglement present in the system. In this work, we consider two types of noise: damping of the oscillator motion in terms of phonon decay, which we denoteκ = κ/ω m , and the number of thermal phonons N th present in the system.
To determine the effect of noise on the entanglement, we evolve the system through equation (14). As expected, we find that the entanglement decreases with time τ as the systems decoheres. We plot the effects of noise on E N in figure 3. In figure 3(a), we have plotted E N as a function of time τ for a noisy environment for differentκ atα = −0.4, z = 1 and N th = 0. Similarly, in figure 3(b), we have plotted E N as a function of rescaled time τ for different squeezing values z atα = −0.4 and N th = 0. We note that while increasing the squeezing z causes E N to increase at first, higher squeezing rates also makes the system more sensitive to noise, a fact that finds extensive confirmation in the existing literature [33][34][35].

Steady state entanglement
After long times τ 1 the system enters a steady state where the entanglement remains at a fixed value. We can obtain this state by solving the matrix equation Aσ (∞) + σ (∞) A T + D = 0 to find the steady state σ (∞) , which is defined by the property All other elements follow by symmetry from the fact that (σ (∞) ) T = σ (∞) , and σ (∞) 11 = σ (∞) 33 . Furthermore, σ (∞) 22 = σ (∞) 44 , due to the symmetry between the two systems. We find the following quantities from which it follows that where Λ = (4α 2 + 8α +κ 2 + 4)/(8α +κ 2 + 4). We note that the number of phonons N th in the system has a strongly detrimental effect on the entanglement. For E N to be maximal, we require that ν − is small. However, sinceν − ∝ N th , the entanglement decreases as N th increases. The decoherence rateκ, on the other hand, is not as influential as the phonon number. Asκ → ∞, we find thatν − = |1 + 2N th |, while as N th → ∞, we find ν (∞) − → ∞. We therefore conclude that reducing the number of phonons in the system takes priority over reducing the decoherence rate. Consequently, the product N thκ which is often quoted in experimental contexts is not as enlightening as an indicator of overall noise levels here, as the two quantities contribute differently to the entanglement.
We further explore the entanglement in the steady state by plotting the logarithmic negativity E (∞) N for noisy dynamics in figure 4(a) for a range of interaction strengthsα ∈ (−0.4, 0.4). As expected, whenκ increases, the logarithmic negativity E (∞) N decreases. The stronger the coupling, the more resilient the system is to noise. Furthermore, all elements in equation (15) are proportional to the term (2N th + 1). To see how this affects the entanglement, we return to the PPT criterion. In figure 4 N as a function of N th andκ. The low values of N th shown on the x-axis give an indication of how sensitive the system is to the number of thermal phonons in the system.

Discussion
We have shown that entanglement from a central potential can be modelled to leading order with Gaussian states in the continuous variable framework. It remains to determine whether such entanglement can be detected in the laboratory.

Error propagation
To what precision must the entries in the covariance matrix σ be determined in order to verify the detection of entanglement? To estimate the required precision, we consider the asymptotic entanglement of the steady state σ (∞) . We then assume that each element σ (∞) i j in the covariance matrix σ (∞) can be determined to within a specific precision i j , which we take to be a percentage of the total value of σ (∞) i j , such that σ (∞) 11 → σ (∞) 11 (1 + 11 ). We also assume that the errors respect the symmetry of σ (so that, for example, 21 = 12 ). From this assumption, we compute the error δE (∞) N via the standard error propagation formula: where in our case k ∈ (11,12,13,14,22,23,24,33,34,44).
To compute the precision required for entanglement detection, we make the sweeping assumption that all errors occur with the same magnitude, such that i j ≡ for all i, j. We then plot the relative error ΔE N = δE (∞) N /E (∞) N in figure 5. From the plots, it becomes evident that detecting logarithmic negativity from small couplings |α| 1 requires a very small percentage error .

Entanglement through central potentials
We now specialise to the following potentials: the Coulomb potential and the Newtonian potential. Given our definition of α in equation (10), the couplings become, respectively: where q 1 and q 2 are the charges on the bipartite system (where the signs will determine whether the potential is attractive or repulsive), 0 is the vacuum permittivity, and G is Newton's constant.
In what follows, we will consider state-of-the-art parameters for each case and determine whether entanglement due to this interaction can be realistically detected. This amounts to examining how close we can setα ∼ −0.5 for each interaction, which is the largest allowed value for an attractive potential.
where e is the electron charge, which results in an attractive potential withα < 0. The number of charges on each sphere can be controlled to exquisite precision by using ultraviolet light [36]. Linearised Coulomb interactions are commonly considered in trapped ions, and have also led to the generation of entanglement [37] but have not yet been implemented for optomechanical systems, although a protocol to enhance cavity-mechanical entanglement through additional Coulomb interactions was proposed in [38]. With the parameters in table 1, it is possible to achieve a couplingα Cl ≈ −0.231, meaning entanglement due to direct Coulombic interactions should be readily implementable for levitated systems. For a different choice of parameters, the attractive coupling can be made even stronger. We note that the mass-scaling is inverted, such thatα Cl ∝ m −1 , which means that entanglement due to the Coulomb potential is suppressed for systems with a larger mass. The number of thermal phonons in the system must still be low at N th = 10 −2 to allow for detection of entanglement, but the decoherence damping rate can be large atκ = 1 as it does not affect the system as much. If we require the error of the entanglement to be no more than 7%, we require that the quantities be determined to within 1% of their value.

Newtonian potential.
We first computeα for the parameters suggested in [19] to see whether gravitational entanglement can be detected with only Gaussian states. With m = 10 −14 kg, r = 200 × 10 −6 m and ω m = 10 Hz, we find α Nw = −8.34 × 10 −16 . Such a weak coupling will not yield any detectable entanglement within this scheme. Even if these figures were more lenient, there is the added complication that the Casimir-Polder interaction will typically dominate over gravity for most parameter choices. Nevertheless, we computẽ α Nw from the optimistic parameters found in table 2. For these values, we still only findα Nw = −6.67 × 10 −8 . For entanglement to be detectable at this interaction strength, we require that as few as 10 −9 thermal phonons are in the system, but we find a slightly more forgivingκ = 0.1. Finally, we require an extremely high precision of = 10 −8 (which we recall is the percentage of the covariance matrix elements). For these values, we obtain an extremely small logarithmic negativity of E (∞) N = (3.3 ± 0.7) × 10 −8 , where the error stands at 22%. Given these numbers, we conclude that gravitational entanglement is not likely to be detectable in the near term with the use of only Gaussian resources.

Experimental detection
A crucial step in the experimental generation of entanglement is the detection stage, where the state is measured to verify the entanglement. Such a detection scheme has, for example, been proposed in the context of a pulsed optomechanical setup [39]. Here, the light-matter interaction is confined to an optical pulse of a timescale much shorter than that of the mechanical oscillation period [40], which means that it can be treated as a single unitary operator that entangles the light and mechanics. By then performing the appropriate measurements on the light, the mechanical state can be inferred by the pulse. In [39], a measurement scheme is proposed to measure all first and second moments of two mechanical elements, which have been previously entangled through optical pulses correlated via the inclusion of a beamsplitter. The entangling step can be readily replaced by the central-potential entanglement scheme considered here.

Conclusions
In this work, we computed the leading-order entanglement due to a generic central-potential interaction between two levitated nanobeads. We derived the Hamiltonian matrix for a linearised potential and investigated the entanglement arising between two initially squeezed states given unitary and noisy dynamics. Furthermore, we derived an analytic expression for the steady state of the system in the noisy setting and proposed a simple continuous-variable test for detecting the entanglement. With these tools, we computed the entanglement from an attractive Coulomb potential and the Newtonian potential. By considering errors that occur when determining the covariance matrix elements, we determined the measurement precision required in each scenario.
Most importantly, we emphasise that this particular setup will not suffice for the detection of entanglement due to gravity. Our results suggest that the inclusion of non-Gaussian resources may play a significant role here (they may, for example, explain the viability of the scheme with freely-falling systems envisioned in reference [19], which relies on the creation of highly non-Gaussian initial states). To prepare two trapped mesoscopic systems in non-Gaussian states, one may utilise the nonlinear optomechanical interaction [41], which couples the photon number to the position of the mechanical element through the addition of a cavity [42]. Specific schemes for generating mechanical cat-states through this interaction have already been proposed by means of the nonlinear optomechanical evolution [43,44], a simple optical interferometry setup [45] or within the pulsed optomechanical regime [46]. We further limited our investigation to Markovian environmental noise in this work, however non-Markovian noisesources are expected to affect optomechanical resonators at low temperatures [47]. It is therefore imperative to investigate the effect of non-Markovian noise on entanglement in this setting, which can be done with a non-Markovian master equation [48]. Finally, to mitigate the small entanglement rates, one may also consider the addition of feedback techniques, which by themselves cannot generate entanglement, but which can serve to amplify already present interactions. These questions, and further scenarios that take into account the noise from the trapping laser will be considered in future work.
After this work was first completed, the authors became aware of similar work by Krisanda et al [49].