Hydrogenic entanglement

Is there any entanglement in the simplest ubiquitous bound system? We study the solutions to the time-independent Schr\"odinger equation for a Hydrogenic system and devise two entanglement tests for free and localised states. For free Hydrogenic systems, we compute the Schmidt basis diagonalisation for general energy eigenstates, and for a Hydrogenic system localised to a three-dimensional Gaussian wavepacket, we demonstrate that measuring its second moments is sufficient for detecting entanglement. Our results apply to any system that exhibits Hydrogenic structure.


I. INTRODUCTION
Are the electron and proton in a Hydrogen atom entangled? In 1926, Erwin Schrödinger successfully predicted the spectral energies of Hydrogen by solving the waveequation for an electron wavefunction in a potential-well created by the positively charged nucleus [1,2]. While the Hydrogenic solutions to Schrödinger's equation have now been known for almost a century, the question of whether the two subsystems are entangled has hardly been investigated. This is likely due to the prevalence of the Born-Oppenheimer approximation [3], which explicitly assumes that the motion of the two subsystems can be treated separately. As a result, the global system factorises into two separable states, and the entanglement is, by definition, removed.
In this work, we go beyond this approximation in order to study the entanglement of the exact Hydrogenic solutions. The ubiquity of Hydrogen in physical, chemical, and biological systems makes it is one of the most well-studied physical systems. Irrespective of the question as to whether entanglement -if present in Hydrogenic systems, has any applications -it is important to know whether such a basic textbook entity of physics is intrinsically entangled.
In recent years, the study of quantum information processing tasks has demonstrated the importance of entanglement to quantum computing [4], quantum cryptography [5], and quantum sensing [6]. A number of fundamental questions, such as whether Nature is described by collapse theories [7], or whether gravity is a quantum force [8], are aided by the quantification and detection of entanglement [9][10][11][12][13]. The ability to experimentally determine entanglement for Hydrogenic systems could have implications for the application of Hydrogenic atoms, artificial atoms [14,15], exiton states [16,17], or possibly even mesoscopic systems, such as levitated nanospheres [18], potentially interacting via a central potential [19].
A rudimentary step in certifying the entanglement in Hydrogenic systems has been taken by Tommasini et al. [20], where they showed that entanglement in a free Hydrogenic ground state system can be verified by diagonalising the state through the Fourier transform and computing its Schmidt coefficients. In the laboratory, however, the assumption that the Hydrogenic system is free is no longer accurate. Rather, any prepared state will be localised to a finite spatial volume, and thus can no longer be diagonalised by the Fourier transform. We must therefore devise an appropriate entanglement test that holds even for localised states.
In this work, we perform the Schmidt basis diagonalisation for arbitrary energy eigenstates of the Hydrogenic system with a free centre-of-mass wavefunction. We find that the spread in the Schmidt basis scales with n −1 , where n is the principal quantum number, which implies that systems in highly excited states are less entangled. These results apply to any system that can be approximated as free, which includes those confined to shallow traps with large centre-of-mass wavefunctions. To treat localised systems, we assume that the centreof-mass wavefunction is prepared in a three-dimensional Gaussian wavepacket. We then take inspiration from the study of entanglement in Gaussian continuous variable states [21,22] to provide a sufficient entanglement test for localised systems. Our results are valid for any energy eigenstates, and succeeds in detecting entanglement for a large number of parameter configurations.
This work is structured as follows. We begin by introducing the Hydrogenic solutions to the time-independent Schrödinger equation in Section II. We then present the Schmidt basis diagonalisation for arbitrary energy eigenstates in Section III for a free Hydrogenic system, and then proceed in Section IV to define a sufficient entanglement test for a Hydrogenic system localised to a three-dimensional Gaussian wavepacket. The work is concluded by a discussion in Section V and some concluding remarks in Section VI.

II. THE HYDROGENIC SOLUTIONS
We begin by considering two systems with a joint wavefunction Ψ(r 1 , r 2 ) that is parametrised by the position vectors r 1 and r 2 . We allow the two systems to interact via a central potential of the form V (r) = α/|r|, where r = r 1 − r 2 and α > 0 is a generic coupling constant that depends on the interaction 1 . The time-independent Schrödinger equation that describes the system is given by (1) As is well known, this equation can be split into a relative and a centre-of-mass equation resulting from the following change of coordinates: where m 1 and m 2 are the masses of the particles. Note that the centre-of-mass and relative momenta P = (p 1 + p 2 ) and p = (m 2 p 1 − m 1 p 2 )/(m 1 + m 2 ) are conjugate to R and r. See Appendix A for details on how the separation of variables is performed. The separated Schrödinger 1 We only consider attractive potentials, which means that V (r 1 , r 2 ) is positive when the minus sign is explicitly included in the Schrödinger equation in Eq. (A2) equations become where µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass of the two-particle system, and where M = m 1 +m 2 is the total mass. The two decoupled time-ndependent Schrödinger equations allow for solutions with a relative wavefunction ψ nlm (r), where n, l, and m are the principal, angular and magnetic quantum numbers, and a centre-of-mass wavefunction ϕ(R) that is an infinite plain wave. Later on, we shall also consider the more realistic case of a localised wavepacket for the centre of mass degree of freedom, and study how the entanglement will be affected by its finite width.
The full state of the system can now be written in terms of the relative coordinate r and centre-of-mass coordinate R as the following separable state: The Hydrogenic wavefunction ψ nlm (r) is given by where L nl (r) is a function of the Laguerre polynomials, and Y m l (θ, φ) are the spherical harmonics. These functions are described in detail in Appendix A.
While the state in Eq. (5) is separable in the {r, R} basis, it cannot be written as a separable state in the original {r 1 , r 2 } basis. As a result, the two subsystems are entangled, and it remains to determine how such entanglement can be detected.
Before we proceed to examine the entanglement of the system, we present the momentum representation of the Hydrogenic wavefunctions, since these will be of great use to us later. They correspond to the Fourier transform of ψ nlm (r), and are given by [23]: where F nl (k) is a function of the Gegenbauer polynomials [24], and Y m l (θ, φ) are again the spherical harmonics, and (k, θ, φ) are the wavevector coordinates. More details on this solution and especially the Gegenbauer polynomials can again be found in Appendix A.

III. ENTANGLEMENT OF THE FREE HYDROGEN ATOM
We start by considering the case when the centre-ofmass wavefunction ϕ(R) in Eq. (5) corresponds to that of a free wavefunction (see Figure 1a): where V is the free space normalisation volume, which we take to infinity at the end of each calculation, R is given in Eq. (2) and K = (k 1 + k 2 ) is the centre-of-mass wavevector of the state. The full eigenstate, including the Hydrogenic part is then given by For all K, one can perform a Galillean transformation into a frame where the centre-of-mass is stationary, corresponding to the product of local unitary operations e im1K·r1/(m1+m2) ⊗ e im2K·r2/(m1+m2) . As one should expect, the entanglement cannot therefore depend on K, which can be set to zero without loss of generality. This leaves us with the following state: In the remainder of this section, we compute the Schmidt basis associated with these eigenstates and generalise the results of Ref. [20] to arbitrary energy eigenstates. We begin by setting out a few preliminaries concerning the Schmidt basis decomposition for continuous variable systems.

A. The Schmidt basis
A necessary and sufficient test to verify whether a pure state is entangled or not consists in casting it in its Schmidt basis and check whether its Schmidt rank is greater than one (in which case the state is entangled). The Schmidt basis and the Schmidt rank are easy to introduce for discrete systems, but have also been studied for continuous variable systems in full generality [25]. We here provide a short review of the Schmidt diagonalisation for discrete systems and show how it can be generalised to continuous systems, with particular reference to the wavefunction notation we are adopting.
It may be shown that a choice of local bases exist such that an arbitrary bipartite state |Ψ can be expanded in a basis formed by pairs of distinct, orthonormal local vectors, which we shall call |w k = |u k ⊗ |v k : where the Schmidt coefficients c k may be taken to be positive and real. The local Schmidt bases are nothing but the eigenbases of the local density operators. It is interesting to note that the local state of a subsystem, defined asˆ 1 = Tr 2 [|Ψ Ψ|] = k c 2 k |u k u k |, satisfies the following eigenvalue equation: Sometimes it is easier to solve this equation for the basis |w k , and thereby reconstruct the state in Eq. (11) than to perform the usual Gram-Schmidt diagonalisation procedure. The link between wavefunciton notation and bra-ket notation is the inner product, where the wavefunctions are scalar functions of the state space, allowing to represent the latter in the (improper) eigenvectors of the position operator. For example: where |x are the position eigenstates and |ψ is the state of the system. In order to introduce the continuous analogue of the Schmidt basis note that, given the wavefunction for a bipartite state Ψ(x, y), where x and y denote the positions of the subsystems (assumed here to be one-dimensional, for simplicity), the analogue of the density matrix is given by which is normalised as follows: dx dy (x, y, x, y) = 1.
The density matrix in the position representation 1 (x , x) of the traced-out subsystem is then given by We concluded before that finding the Schmidt basis for discrete states amounts to solving the eigenvalue equation in Eq. (12). The continuous analogue of the eigenvalue problem is the following integral equation, otherwise known as a Fredholm equation of the second kind [26]: where {λ i } is a set of eigenvalues ofˆ 1 parametrised by the index i. Since we in this work are interested in bipartite systems in three spatial dimensions, we write the density matrix of the full state as This allows us to rewrite the Fredholm equation in Eq. (17) in terms of a three-dimensional subsystem state 1 (r 1 , r 1 ), given by whence the Fredholm equation becomes where we have indexed φ k with the (possibly) continuous index k. Furthermore, in terms of the full state, the Fredholm equation in Eq. (17) becomes dr 1 dr 2 Ψ * (r 1 , r 2 )Ψ(r 1 , r 2 ) φ k (r 1 ) = λ k φ k (r 1 ) . (21) To obtain the diagonalisation in three dimensions, Eq. (20) must be solved for φ k . We anticipate that finding the solution will in general be an extremely challenging task. In the next Section, however, we show that Eq. (20) has a simple solution when the system is invariant under spatial translations.

B. Schmidt basis diagonalisation of a translationally invariant state
Before considering the Schmidt diagonalisation of the free Hydrogenic wavefunctions, we proceed with a proof showing that translationally invariant ("homogeneous") states are diagonalised by the Fourier transform, which was first demonstrated in Ref. [20]. We define translationally invariant by the fact that displacing the two subsystems by the same distance in space leaves the wavefunction invariant, up to a global phase.
A translationally invariant wavefunction must depend on the variable r alone, and not on R, since the former is invariant under an equal translation of r 1 and r 2 whilst the latter may be modified arbitrarily by such translations.
This implies that the local density function 1 (r 1 , r 1 ), whose eigenbasis determines the Schmidt decomposition, is a function only of the difference of its arguments, since Eq. (19) in this case reads = dr 2 Ψ * (r 2 )Ψ(r 1 − r 1 + r 2 ) = 1 (r 1 − r 1 ), where we have made the variable substitution r 2 → r 1 − r 2 , treating r 1 as a constant. The Fourier transform and inverse of 1 (r 1 − r 1 ) may then be expressed in terms of a single (three-dimensional) variablẽ where we have defined r = (r 1 − r 1 ) and its conjugate variable k = (k 1 − k 1 ). 2 We then recall the Fredholm eigenvalue equation for continuous systems given in Eq. (20). Inserting the translationally invariant state and the momentum eigenstate ansatz φ k = e ik·r 1 , we find where we have used the Fourier transform of (r) shown in Eq. (24). This demonstrates that the local density operator is diagonal in the local momentum eigenbasis, which thus form its associated Schmidt basis, with Schmidt coefficients given by the Fourier transform˜ (k): a translationally invariant state is always diagonalised by the Fourier transform.

C. Entanglement for general Hydrogenic energy eigenstates
We have shown that the local density operator of homogeneous state is always diagonalised by the Fourier transform. A free Hydrogenic system, whose wavefunction is shown in Eq. (10), fulfills the criterion of translation invariance since it depends solely on the distance r = (r 1 − r 2 ) between the two systems.
In order to compute the Schmidt coefficients of a Hydrogenic eigenfunction, and hence qualify its entanglement, we need first to trace out one of the two particles, as follows 1 (r 1 − r 1 ) = dr 2 Ψ * (r 1 , r 2 )Ψ(r 1 , r 2 ) Here, through the substitution y = r 1 −r 2 , we have again shown that the traced out Hydrogenic state is translation invariant, since it only depends on the local variable r = r 1 −r 1 . As we proved in the previous Section, the Fourier transform˜ 1 (k) corresponds to the Schmidt basis of the a Hydrogenic system, where k = k 1 − k 1 is the Fourier complement to r. The state is separable if and only if the probability distribution˜ (k) is a delta-function, which provides one with a criterion to test entanglement in this context. This corresponds to the discrete analogue of having one single non-zero Schmidt coefficient.
Notice that, even if the continuous analogue of the Schmidt coefficients are known, an actual quantification of the entanglement is difficult in this case. The von Neumann entropy is in fact not a good quantifier for such continuous systems, as it may be shown that a state with diverging von Neumann entropy may be found, in the trace norm topology, arbitrarily close to any state (nor does the finite value of the free hydrogenic Hamiltonian heal this state of affairs, since its spectrum is continuous) [27]. The local linear entropy S Lin = 1 − Tr[ˆ 2 1 ] does however converge, but its evaluation for general eigenstates is unwieldy and does not admit a direct operational interpretation. We return to the linear entropy at the end of this section.
We opt instead to evaluate a heuristic quantifier of entanglement, already suggested in [20], given by the standard deviation of the Schmidt function: This does in a sense quantify the deviation from the delta function that characterises separable states: if ∆k = 0, the state is separable, and if ∆k > 0, the state is entangled. We proceed to compute ∆k. The continuous analogue of computing expectation values with the trace operation is given by We then insert the Fourier transform in Eq. (24) to find where V appears because we integrate over all space. This expectation value k is calculated by integrating over all vectors, which averages to zero: The variance k 2 , on the other hand, is given by Using the momentum representation of the Hydrogenic wavefunctions, shown in Eq. (7), and realising that the free system is proportional to V −1 , we find The integrals over θ and φ satisfy the normalisation condition for the spherical harmonics (see Eq. (A19) in Appendix A), and we are left with The integral overk returns a well-known result from atomic physics. It admits the standard solution [23]: We prove this relation, which extends what was previously known for the Hydrogenic ground state to any energy eigenstate [20], in Section B 2 in Appendix B.
We conclude that the standard deviation of the local momentum, related to the local wavevector byp = k , reads where the expression for the reduced Bohr radius a 0 = 2 /(αµ) was inserted (recall that µ is the reduced mass of the system and that α is the potential interaction strength, which has dimensions of an energy times a length).
The standard deviation of the continuous local spectrum only depends on the principal quantum number n, and not on l and m. This is the case since the operations that transform sectors with the same n and different l and m are local unitary operations (which is particularly clear for the azimuthal number m, whose different values are related by rotations), and thus cannot affect properties related to the Schmidt spectrum and entanglement.
On the other hand, the distribution of local momentum is more spread out for lower principal quantum numbers, which is non-trivial and suggests that, in practice, lower quantum states would be more favourable for the detection of quantum correlations in Hydrogenic systems. Notice that the fact that we can qualify the presence or absence of entanglement so liberally (by any positive value of the standard deviation!) is just an artefact of the fact that our theoretical finding applies to ideal (pure) eigenstates.
Further, our quantifier grows with the interaction strength, as one should expect, and with the reduced mass, which indicates that, at a specific given total mass, the entanglement is maximum for two equally distributed masses. However, the eigenstate of a Hydrogen atom would be about twice more entangled than the corresponding positronium eigenstate, since the large mass of the proton nucleus factors out of µ and leaves µ = m e , where m e is the electron mass. Compare this with the reduced mass of the positronium eigenstate (a bound state of one electron and one positron), which reads µ = m e m p /(m e + m p ) = m e /2. We also note that the expression in Eq. (35) is not independent of the choice of Fourier normalisation, however, as this is merely a scaling factor, the choice does not significantly affect the results.
Before we proceed to investigate localised Hydrogenic states, let us briefly return to the linear entropy S Lin of a Hydrogenic system. As mentioned above, the linear entropy converges for a free Hydrogenic systems, although it does not admit a direct operational interpretation as a measure of the entanglement. We provide this result for completeness.
In Appendix C, we show that the linear entropy S Lin for a free Hydrogenic system is given by While we find that the integral in Eq. (36) converges, we conclude that when V → ∞, the linear entropy tends to S Lin = 1, which conventionally indicates that the state is maximally entangled. However, as mentioned before, we primarily include this result for completeness, because a direct operational interpretation of the state as entangled as a result is not necessarily valid.

IV. ENTANGLEMENT OF A LOCALISED HYDROGENIC SYSTEM
The approximation of the system as free breaks down when the overall centre-of-mass wavepacket is of a similar width compared with the characteristic length scale a 0 of the system. For this setting, the methods and results demonstrated in the previous section are no longer valid if the system is prepared in a localised state in the laboratory.
By definition, a localised state is no longer invariant under spatial translations, which means that it can no longer be diagonalised by the Fourier transform. Finding the Schmidt basis of the state would require solving the integral equation shown in Eq. (20), which we anticipate to be an extremely difficult task.
In this work, we choose instead to explore an entanglement test based on detecting a violation of the positivity of the partial transpose through the information contained in the second statistical moments of the localised Hydrogenic state [21]. One should bear in mind that such a test applied to non-Gaussian states is only sufficient, and not necessary for entanglement.
We consider therefore the case where the centreof-mass wavepacket a three-dimensional Guassian wavepacket centred at the origin, as shown in Figure 1b. While a Gaussian wavepacket is not a solution to the time-independent Schrödinger equation of a free particle, it is nonetheless a very reasonable localised state, which can be prepared and maintained, for instance, by trapping the system in a harmonic potential that addresses the centre-of-mass degree of freedom. The addition of dynamics into this setting is a problem we leave for future work.
The Gaussian wavefunction is given by where b has units of length and encodes the standard deviation of the wavepacket at the specific moment in time for which we consider the system. The full state of a localised Hydrogenic system therefore becomes It is clear that the state appears entangled even in the Gaussian wavepacket alone, since ϕ(R) in Eq. (37) cannot be written as a product of functions of r 1 and r 2 . We therefore expect contributions to entanglement from both the relative and centre-of-mass degrees of freedom.
We proceed now to study the entanglement criterion for this state. For simplicity, we restrict the analysis to the case of equal masses, such that m 1 = m 2 .

A. Entanglement in the first and second moments
In order to devise an entanglement test for a localised Hydrogenic state, we can adopt well-established techniques developed in the study of quantum information with continuous variables [28]. In particular, a necessary and sufficient entanglement criterion based on the second order statistical moments of the canonical operator exist for two-mode Gaussian states [21,29], which is however sufficient to detect the entanglement of any quantum state. 3 Such a criterion is based on a general necessary test for the positivity of the partial transpose (PPT) at the level of second moments, and on the observation that, in turn, PPT is necessary for the separability of arbitrary states. A necessary criterion for separability is equivalent to a sufficient criterion for its converse, namely quantum entanglement.
An entanglement test on the second moments can be carried out as follows. Given a suitable canonical basis of self-adjoint operatorsr = (x 1 ,p 1 ,x 2 ,p 2 . . .x n ,p n ) T for n modes, the covariance matrix σ is given by whereˆ is a Gaussian (or non-Gaussian) state, and where the transpose T is taken in the outer-product sense, including all possible pairs of operators to form a real, symmetric matrix of correlations. For continuous variable systems, partial transposition is equivalent to changing the sign of the second canonical variable. This implies that half of the offdiagonal elements in the covariance matrix obtain a minus sign [28]. Because any physical (i.e., derived from a trace-class, positive density operator) covariance matrix σ must satisfy the uncertainty principle that can be cast as −(Ωσ Tp ) 2 ≥ 1, where Ω is the symplectic form defined in this basis as the sufficient entanglement test can be stated as whereν j are the N symplectic eigenvalues of the partially transposed covariance matrix σ Tp . The symplectic eigenvalues are defined as the square roots of the eigenvalues of −(Ωσ Tp ) 2 , which are at least two-fold degenerate (so that there are N independent ones in the case on hand). Hence, in order to demonstrate entanglement for the Hydrogenic wave-function (which are globally non-Gaussian due to the relative coordinate contribution, even if the centre-of-mass coordinate is taken in a Gaussian wave-packet, which would amount to a Gaussian state in the quantum optics nomenclature), it suffices to show that one of the symplectic eigenvalues of its covariance matrix satisfiesν j < 1. We therefore need to evaluate its second moments, which is what we do in the following section.

B. Second moment entanglement for Hydrogenic states
The bipartite state in Eq. (38) has support in three spatial dimensions. This setting yields six independent modes (one for each spatial direction for each of the subsystems). As a result, the covariance matrix is a 12 × 12 matrix (with two variables per mode). To construct the matrix, we must compute the expectation values and variances for the position and momentum variables in each mode.
In the original partition in terms of the r 1 and r 2 coordinates, we construct the following (dimensionfull) basis of coordinates: which partitions the covariance matrix in terms of the two original subsystems.
It would be cumbersome to compute the first and second moments directly in the basis X of Eq. (42), since the variables r 1 and r 2 are mixed between the Hydrogenic eigenstates and the Gaussian wavepacket. However, we can simply evaluate the statistical moments in the {r, p, R, P} basis, separately for the Hydrogenic wavefunction and the Gaussian wavepacket, and then apply the symplectic transformation S that relates the local to the decoupled coordinates on the covariance matrix itself, which transforms, by congruence, as σ = SσS T . With the S shown in Eq.(D3) in Appendix D, the transformed basis vector in Eq. (42) becomes where x, y, z and p x , p y , p z are the relative position and momentum coordinates given by (for equal masses) and so on. Similarly, X, Y, Z and P X , P Y , P Z are the centre-of-mass position and momentum coordinates given by and so on.
We are now ready to compute the first and second moments of the Hydrogenic states. To ensure that they are dimensionless, we rescale the position coordinates by the characteristic length scale, the reduced Bohr radius a 0 (see Appendix A), and the momentum coordinates by the characteristic momentum /a 0 .
Using expressions for the Legendre polynomials and known relations such as the Kramer-Pasternak relation [30,31], we obtain the following first and second moments in the {r, R} basis. See Appendix E for the full, rather lengthy calculation. The dimensionless expectation values and variances for the relative position variables of the hydrogenic eigenstates are and, since x 2 /a 2 0 = y 2 /a 2 0 : The expectation values of the relative momentum variables are and the variances are given by, with p 2 x a 2 0 / 2 = p 2 y a 2 0 / 2 : As for the Gaussian wavepacket, the centre-of-mass expectation values and variances read We now promptly transform back to the original basis X = S T X , and σ = S T σ S. We then apply the PPT criterion and compute the symplectic eigenvaluesν j of σ Tp . They are given by, in terms of the previously computed variances, Due to the rotational symmetry of the system, the symplectic eigenvalues are degenerate withν 1 =ν 3 =ν 5 and thatν 2 =ν 4 =ν 6 . The two unique eigenvalues are given byν To detect entanglement, we merely require that one of the eigenvalues becomes smaller than unity: ν (j) nlm < 1 for j = 1, 2. The first symplectic eigenvalueν (1) nlm indicates the entanglement increases with b. In fact, as b → ∞, ν j → 0, which implies that the logarithmic negativity for the second moments diverges for a free Hydrogenic system.
The second eigenvalueν (2) nlm has the opposite dependence of a 0 and b, but scales as n −1 , which means that the larger n becomes, the more entangled the state appears. We also note that both eigenvalues depend on the ratio a 0 /b, and that the dependence of the magnetic quantum number m is purely detrimental to entanglement detection; the larger m is, the less likely we are to detect the entanglement.
For the Hydrogenic ground state, with n = 1 and l = m = 0, the eigenvalues becomẽ These values depend purely on the ratio of the Hydrogenic characteristic length-scale and the Gaussian wavepacket spread a 0 /b. However, entanglement can be detected for all values through examination of these two eigenvalues. When a 0 ≤ b, the second eigenvalue always detects entanglement. For the special case when b/a 0 = 3/2, entanglement can briefly not be detected at all, until b/a 0 ≥ 5/3. We have plotted the function min(ν 100 ,ν 100 ) in Figure 2 to demonstrate the regimes where entanglement can be detected. The blue area indicates where either of the eigenvalues dip below 1, and where entanglement can be inferred. The light yellow area indicates where both eigenvalues are larger than 1, where entanglement cannot be inferred.
We conclude that in order to apply this entanglement test, one must measure the relative and centre-of-mass coordinates in a single spatial direction. That is, it suffices to measure x 2 , p 2 x , X 2 and P 2 X to compute the two symplectic eigenvalues in Eq. (52). It is worth mentioning that this test, in contrast to the the previous one based on the relative momentum variance, would be sufficient also for noisy, mixed states (albeit for such states it would likely be more difficult to violate).

V. DISCUSSION
We derived two expressions that can serve as entanglement tests for a free Hydrogenic system (see Eq. (35)) and a localised Hydrogenic system (see Eq. (52)), respectively. Here, we briefly discuss our results and their potential application to detecting entanglement between two subsystems that interact via a central potential.
The first additional system that might exhibit Hydrogenic structure is positronium, which is the bound state of one electron and one positron. With a half-life of 0.1244 ns [32], it is a highly unstable system. We note that detecting entanglement in positronium might be extremely challenging, however it may also constitute one of the most straight-forward routes towards detecting entanglement between matter and anti-matter systems.
We speculate that another candidate family of systems that might exhibit Hydrogenic structure could be mesoscopic systems that interact through a central potential, such as a Coulombic potential. Recent advances in the ability to control the number of charges on levitated silica spheres [33] or implant changes through the addition of nitrogen-vacancy (NV) centres [34] provide an excellent means of control. It is currently an open question as to whether mesoscopic systems would at all readily exhibit Hydrogenic structure, but should this be the case, the methods developed here may be used towards establishing any entanglement present between them.
We further note that any atom with Hydrogenic structure, such as the Helium ion, or any two oppositely charged interacting ions could be treated with our results. Artificial atoms in semiconductor structures are yet another possible setting where tests similar to the ones proposed here might become of relevance. In exitonic bound states, one can plausibly have access to each of the entangled components to examine their entanglement.
Finally, one may also look for entanglement in phenomena beyond the Standard Model. These include proposals for the existence of so-called magnetic monopoles [35], which, like opposite electric charges, are theorised to form bound systems of north-south magnetic monopoles. These are referred to as monopolium [36]. Should these particles be observed, their interaction might cause the two monopole subsystems to become entangled in a similar fashion to the Hydrogen atom.

VI. CONCLUSIONS
In this work, we proposed two entanglement tests for free and localised bipartite systems with Hydrogenic structure. We computed the spread of the Schmidt spectrum for an arbitrary energy eigenstate of a delocalised (free) Hydrogenic system, and showed that it scales with the inverse of the principal quantum number n −1 . For localised systems, we showed that measuring the variances of the relative and centre-of-mass coordinates in a single spatial direction suffices to demonstrate that the system is entangled. We believe that these results could potentially aid the study of entanglement for systems that display Hydrogenic characteristics, such as entanglement between matter and anti-matter systems or freely-falling mesoscopic systems.
A number of open questions remain, such as the extension of these results to dynamical processes, including the formation of Hydrogenic states starting from individually localised systems. We leave these questions to future work.
In this Appendix, we introduce the Hydrogenic solutions to the time-independent Schrödinger equation in three dimensions for two particles that interact through a central potential. We describe here the general procedure for separating the Schrödinger equation in terms of the centre-of-mass R and relative coordinate r, which will be used to describe the interaction. We begin by presenting the procedure in position space, and then move on to review the corresponding solutions in momentum space.

Hydrogenic wavefunctions in position space
We follow the standard derivation presented in Appendix 8 of Ref. [23]. The Schrödinger equation for N interacting particles in first quantisation with coordinates r i and masses m i , is given by where m j is the mass of particle j, ∇ 2 j is the Laplacian operator, and where V (r 1 , r 2 , . . . , r N , t) is the sum of central potentials affecting the particles. For two interacting particles, the equation reduces to where the explicit form of the potential is given by where α > 0 is a coupling constant that depends on the interaction at hand, and r = r 1 − r 2 . For the Coulomb potential, we find that α = q 1 q 2 /(4π 0 ), and for a Newtonian gravitational interaction, α = Gm 1 m 2 . The two-particle Schrödinger equation in Eq. (A2) is solved by moving to a centre-of-mass coordinate R and a relative coordinate r, defined as (with conjugate momenta P = (p 1 + p 2 ) and p = (m 2 p 1 − m 1 p 2 )/(m 1 + m 2 ), defining a canonical quadruple). We then introduce the reduced mass µ = m1m2 m1+m2 and the total mass M = m 1 + m 2 in order to write the original coordinates as It can then be shown that the gradient operators can be written as which allows us to rewrite the two-particle Schrödinger equation as Since the potential only depends on r and not on R, we can separate the variables into the Schrödinger equation above, as per Ψ(r, R) = ψ(r)ϕ(R), obtaining where E CM is the centre-of-mass energy and E rel is the relative energy. The total energy of the system is Eq. (A11) describes the free centre-of-mass wavefunction (the overall wavepacket) of the two-particle system and admits wave-like solutions. Eq. (A12), on the other hand, admits the Hydrogenic solutions ψ(r) with E rel being the usual Hydrogen energy spectrum [37]. The solutions to the relative equation are given by where n, l, m are the principal, angular and magnetic quantum numbers, respectively. r, θ and φ are elements of the relative coordinate vector: r = (r, θ, φ) T , and where L 2l+1 n−l−1 (2r/na 0 ) is a Laguerre polynomial, defined by and where L q (x) is an associated Laguerre polynomial, given by Finally, we note that a 0 is a length scale set by the interaction coefficient α, corresponding to the reduced Bohr radius: Returning to the functions shown in Eq. (A14), the spherical harmonics Y m l (θ, φ) are given by where P lm (cos θ) are the associated Legendre polynomials without the Condon-Shortley phase [26]. The spherical harmonics satisfy the following orthogonality relation We make frequent use of this relation in the following Appendices. This concludes our summary of the Hydrogenic solutions in position space.

Hydrogenic wavefunctions in momentum space
The Hydrogenic wavefunctions in momentum spaceψ(k), where the relation between the wavevector k and the momentum is p = k, are defined as the Fourier transform of the position space wavefunctions: ψ(k) = dr e −ik·r ψ(r). (A20) Note that we denote this k with a bar in the main text to differentiate it from another relative quantity, but we omit this here for notational clarity. The explicit form of the Hydrogenic wavefunctions in momentum space can be computed to be (see Appendix 5 in [23]):ψ where and where C α N (x) are the Gegenbauer polynomials [24], which are defined by the relation with |s| < 1. The Gegenbauer polynomials are special cases of the more general Jacobi polynomials and play a role that is analogous to the role played by Legendre polynomials in the theory of the three-dimensional spherical harmonics [38]. The Gegenbauer polynomials are normalised by the following orthogonality relation where α is a dummy index, not to be mixed up with the coupling constant defined earlier. This is yet another relation that is used many times in what follows. The normalisation of the momentum wavefunctions reads dk |ψ(k)| 2 = 1, which, since the spherical harmonics are orthogonal (see Eq. (A19)), implies that This concludes the summary of the Hydrogenic wavefunctions in momentum space. We return to these expressions in Appendix B, where we compute the variance of the relative momentum variable.

Appendix B: Normalisation and expectation values for the momentum Hydrogenic wavefunction
In this Appendix, we prove Eq. (34) in the main text, which provides a relation for the variance of the relative momentum coordinate. This, in turn, corresponds to the continuous analogue of the Schmidt coefficients of a free Hydrogenic system, and is therefore a key result in this work. We begin by proving the orthogonality relation and the normalisation condition of the Gegenbauer polynomials, as this will become useful in proving Eq. (34).
We used the notation k in the main text to refer to the local wavevector rather than that related to the relative momentum (which is defined as k = (m 2 k 1 − m 1 k 2 )/(m 1 + m 2 )). We here remove the bar for notational convenience.

Normalisation of the Hydrogenic momentum wavefunction
As discussed in Appendix A, the momentum representation of the Hydrogenic wavefunctions is given in Eq. (A21). We begin by proving the orthogonality relation in Eq. (A26), which we reprint here for brevity: We reprint the expression for F nl (k) here for convenience: We note that momentum wavefunctionψ(k) has units of length 3/2 due to the appearance of a 3 0 in Eq. (B2). This is to ensure that the normalisation is dimensionless.
Our goal is to prove the relation in Eq. (B1). We start by writing out the expression in full: We now multiply and divide by (a 0 k) 4 and 4 l+2 in order to collect the following terms: We now define the new variable The variable substitution implies that and the derivative becomes With this substitution, the limits become {−1, 1}, which is what we require for the normalisation condition of the Gegenbauer polynomials. We insert this into Eq. (B4) to find The last integral in Eq. (B8) can now be written as two separate terms by multiplying out the first bracket: The integral in the first term of Eq. (B9) can be evaluated by using the orthogonality relation in Eq. (A24). With the index substitution α → l + 1 and n → n − l − 1, we find (B10) When we include the prefactor, the full first term of Eq. (B9) evaluates to 2 π Since the Gamma function is given by Γ(n) = (n − 1)! for a positive integer n, we are able to write This is the normalisation we require. It remains to show that the term on the left-hand side of Eq. (B11) vanishes. As noted in Ref. [39], to complete the proof, we can use the following recursion formula for the Gegenbauer polynomials: With this relation, the second integral in Eq. (B11) (without its prefactor and with the appropriate index substitutions) can be written as Since the indices of the Gegenbauer polynomials inside and outside the bracket differ, the integral vanishes due to the orthogonality of the Gegenbauer polynomials in Eq. (A24). This concludes our proof of Eq. (B1).

Proof of Hydrogenic wavefunction momentum variance
In Section III C, we showed that the variance of the relative momentum variable k (referred to ask in the text) corresponds to the continuous analogue of the Schmidt coefficients of the system, which is given by the concise expression in Eq. (35). Here, we prove the relation in Eq. (35).
We wish to show that where F nl (k) is given in Eq. (B2). This allows us to write the integral as By now multiplying and dividing by a 4 0 and 4 l+2 , we can write Eq. (B16) as To evaluate the integral, we perform the same variable substitution as in Eq. (B5). Given the relationships in Eqs. (B6) and (B7), and the new limits The last integral in Eq. (B18) can again be written in terms of two separate terms: As in the previous section, the second term vanishes due to the mismatched indices. Then, by using the orthogonality condition for the Gegenbauer polynomials in Eq. (B10), we find that the first term evaluates to (B20) Noting that Γ(n) = (n − 1)!, this expression simplifies to which is the expression we wanted. This concludes the proof of Eq. (35).

Appendix C: Linear entropy of the Hydrogen atom
The goal of this Appendix is to compute the linear entropy of a Hydrogenic system. We provide this result for completeness and we emphasise that it is not possible to provide an operational interpretation of the linear entropy as a measure of the entanglement when the basis is continuous. In contrast, we refrain from deriving an expression for the von Neumann entropy, since apart from the lack of operational interpretation, one cannot compute an entropy of a dimensionfull quantity.
The continuous variable analogue of the linear entropy S Lin is where (r 1 , r 1 ) is the continuous analogue of the subsystem density matrix, given by tracing out one of the variables: (r 1 , r 1 ) = dr 2 Ψ * (r 1 , r 2 )Ψ(r 1 , r 2 ).
The form of the linear entropy in Eq. (C1) is the continuous analogue of the discrete matrix multiplication. Inserting the Fourier transform Eq. (24) of the state (r 1 , r 1 ) into Eq. (C1) gives S( ) = 1 − dr 1 dr 1 (r 1 , r 1 ) (r 1 , r 1 ) where we have defined the Dirac delta function as δ(k − k) = 1 (2π) 3 dx e ix·(k −k) , and With the definition of the momentum space Hydrogenic wavefunctions in Eq. (A21), we find To the authors' knowledge, this integral has no known standard solutions. Instead, we aim to derive a closed-form expression that can be evaluated for specific choices of n, l and m. To simplify the evaluation of the integral, we write it in terms of a radial integral I k and an angular integral I θ , such that We now proceed to evaluate each integral separately.

The angular integral
Let us start by evaluating the angular integral, which involves four spherical harmonics: To evaluate this integral, we first seek to reduce the order of the spherical harmonics from four to two. We use the following Clebsch-Gordon expansion [40] to write where the 3 × 2 matrices are in fact scalar functions known as the Wigner 3-j symbols. They are defined in terms of the Clebsch-Gordon coefficients and read: where j 1 m 1 j 2 m 3 |j 3 (−m 3 ) denotes the inner product between two eigenstates. The Wigner 3-j are zero unless the following conditions are fulfilled: When l 1 = l 2 = l and m 1 = m 2 = m, as in our case, we find We proceed to compute the radial integral.

The radial integral
The radial integral is given by Now, F nl (k) is given in Eq. (B2). Just like we did in Appendix B, we rewrite the integral in Eq. (C16) by dividing and multiplying with a 8 0 k 6 and 4 2l+4 to find: We again perform the same variable substitution as in Eq. (B5) to find:

Final expression for the linear entropy
Combining the angular integral in Eq. (C15) and the radial integral in Eq. (C23), we find This expression is not particularly elegant (perhaps there are also other, more suitable ways to evaluate the radial integral, although the present authors were not able to find a way to do so), but it can be evaluated for specific choices of n, l and m. For the Hydrogenic ground state, for example, the linear entropy becomes where we note that the fraction is dimensionless due to the appearances of V and a 0 . Since we take V → ∞ at the end of each calculation, we are left with S Lin = 1, which conventionally implies that the entropy is infinite (however, as mentioned, the linear entropy cannot be interpreted as a measure of entanglement). A similar result for the ground state linear entropy of the Hydrogen atom was reported in a Masters thesis, see Ref. [42].

Appendix D: Covariance matrix and symplectic eigenvalues
In the main text, we noted that successfully detecting entanglement in the first and second moments of a Gaussian state immediately implies that a non-Gaussian state with the same moments is also entangled [28]. In this Appendix, we demonstrate in detail how we first perform a symplectic transformation from the {r 1 , r 2 }-basis into the relative {r, R} basis. We provide the explicit form of the symplectic transformation and compute the symplectic eigenvalues of the final covariance matrix.
Our task is to compute the entanglement between the two subsystems in the {r 1 , r 2 } basis. We define σ in this basis with elements given by The diagonal elements of σ are given by diagσ = 2 x 2 1 , p 2 x,1 , y 2 1 , p 2 y,1 , z 2 1 , x 2 2 , p 2 x,2 , y 2 2 , p 2 y,2 , z 2 2 , p 2 Where x j , y j , z j and p x,j , p y,j , p z,j are the position and momentum coordinates of the two respective systems.
To transform into the basis of the subsystems parametrised by r 1 and r 2 , we perform a symplectic transformation S that maps the covariance matrix from the {{r, p, R, P}} basis to the {r 1 , p 1 , r 2 , p 2 } basis.
The symplectic transformation S that maps the covariance matrix from the relative and centre-of-mass basis to the subsystem basis is given by In language commonly used in the context of quantum optics, this transformation is equivalent to the combination of a 50:50 beam-splitter and a diagonal squeezer.
The new covariance matrix is given by σ = SσS T . We prove in the following Appendix that all the first moments in the relative basis are zero, and that σ is a diagonal matrix with the following entries: where we have defined the relative and centre-of-mass position coordinates as follows: and, similarly, the momentum variables: p y = p y,1 − p y,2 2 , P Y = p y,1 + p y,2 , Once the expectation values of the coordinates in Eq. (D6) have been computed, we transform back into the original basis with σ = S −1 σ S −1 T . To better see the explicit action of this transformation and the return of coherence to the off-diagonal elements, we write σ in terms of block matrices: Then, since the subsystems are symmetric, with (r 1 , r 1 ) = 2 (r 1 , r 1 ), the covariance matrices for both systems are the same with σ 1 = σ 2 , where σ 1 is given by

21
The off-diagonal elements, which contain the correlations between subsystem 1 and 2, are given by (D9) We are now ready to perform the partial transposition of σ and compute its symplectic eigenvalues to determine whether the r 1 and r 2 subsystems are entangled. We call ν j the symplectic eigenvalues of σ andν j the symplectic eigenvalues of the partially transposed covariance matrix σ Tp .
In a continuous variable system, partial transposition is equivalent to introducing a minus sign to all the momenta (or all the positions) of a subsystem [28]. This implies that every element corresponding to the relative momentum expectation values p x,2 , p y,2 , p z,2 becomes inverted. The diagonal elements of σ Tp remain invariant, and the off-diagonal elementsσ 12 are now given bỹ (D10) We can now compute the symplectic eigenvalues of the partially transposed system. They are given by the eigenvalues of the product iΩσ Tp , and we find: We note that the symplectic eigenvalues do not mix between the modes. Entanglement can then be established once the first and second moments of the Hydrogenic system have been computed, which we do in the next Appendix. Due to the rotational symmetry of the system, these three pairs of quantities will take the same along each spatial direction, so that our entanglement test will effectively reduce to a two-mode problem.
where ψ(r) are the Hydrogenic eigenstates shown in Eq. (A14). We proceed to compute all the elements of the covariance matrix. Given the definitions of the relative and centre-of-mass coordinates in Eq. (D6), full two-mode covariance matrix is a 12 × 12 matrix that is explicitly given by The off-diagonal elements are given by a similar algorithm for i = j. However, since the wave-function is a product of terms that depend separately on each such variable, we only have to compute the diagonal elements. The advantage from working in this basis now becomes evident. Depending on whether the variable of interest belongs to the relative or the centre-of-mass coordinates, we find that the other part of the state does not influence the result. In other words, if we are computing a moment of one of the centre-of-mas variables R j , then we find where R 6 denotes the integration over the full spatial domain of the bipartite state. When the state separates in terms of the relative and centre-of-mass coordinates as Ψ(r, R) = ψ(r)ϕ(R), the part of the state that does not relate to the expectation value in question simply satisfies the normalisation conditions. For example, if we wish to compute the expectation value where the first integral is just the normalisation condition for ψ nlm (r)., we are left with the task to evaluate the Guassian integral, or the integral over the Hydrogenic eigenstates. Finally, we remark that since the transformation S, shown in Eq. (D3), between the {r 1 , r 2 } basis and the {r, R} basis is symplectic, the Jacobian for the substitution of variables in the integrals is equal to unity.

Expectation values and variances for Gaussian centre-of-mass wavepacket
If the Gaussian wavepacket is centred at the origin, its expectation values are zero. This is the case for the state we consider, and thus we find The variances, on the other hand, are given by The first integral is equal to and the second integral is equal to So we are left with This element will be the same for all other quantities. In summary, we have Next, we wish to compute the centre-of-mass momentum expectation values. The momentum variable in the X direction is given by and similarly for P Y and P Z . Taking the expectation value, we find where the integral vanishes because the function is odd. Due to the symmetry of the wavepacket, it follows all the expectation values of the momentum coordinates are zero. We proceed to compute the variance of the momentum variable. We find The first integral again evaluates to that in Eq. (E9). The second integral becomes Putting everything together, we find Again, due to symmetry, we find P 2 X = P 2 Y = P 2 Z . So in summary, the momentum variances become We proceed to consider the Hydrogenic wavefunctions.

Position expectation values of the Hydrogenic subsystems
We wish to compute the expectation values x , y , and z . We recall that the Hydrogenic wavefunctions are given by We can rewrite this as ψ nlm (r) = L nl (r)Y m l (θ, φ). (E20) When we take the absolute value |Y m l (θ, φ)| 2 the result becomes independent of φ. We write this as We proceed to compute the first moments of the position variable. The expression that we must calculate becomes where the cross-terms vanish because of the orthogonality relation in Eq. (E40). Then, we use the same relation to find that In summary, we have that x 2 = a 2 0 n 2 (5n 2 − 3l(l + 1) + 1) 2 l 2 + l + m 2 − 1 (2l − 1)(2l + 3) , y 2 = a 2 0 n 2 (5n 2 − 3l(l + 1) + 1) 2 l 2 + l + m 2 − 1 (2l − 1)(2l + 3) , z 2 = a 2 0 n 2 (5n 2 − 3l(l + 1) + 1) 2 The fact that z 2 differs from the other two variances is reasonable given the identification of the rotation axis.

Variance of px, py and pz
The variances can be quickly computing by realising that we obtain integrals over the spherical harmonics that we already evaluated in Section E 3.
We start with the first variance p x , however we will first compute k x and then multiply by to obtain the momentum. This quantity is easier to compute in the momentum basis. We therefore use the Fourier transform in Eq. (23) and Eq. (24), to write The Hydrogen wavefunction in the momentum representation is given by, in terms of the wavevector k, in Eq. (A21). Using this expression, and moving to spherical coordinates (k, θ, φ), where k x = k cos θ sin φ, we find where we again used the fact that |Y m l (θ, φ)| 2 = Y(θ), which is independent of φ. However, we again note that the last integral is again zero, and thus p x = 0. In fact, since the angular integral is the same for the position and momentum expectation values, we conclude that we also have p y = p z = 0. Intuitively, this is reasonable, since a non-zero momentum expectation value would mean that the system has a net non-zero motion in one of the directions. 6. Variance of p 2 x , p 2 y and p 2 z We proceed to compute the variances of the relative momentum variables. In the wavevector basis, we find that = ∞ 0 dp p 4 |F nl (p)| 2 π 0 dθ sin 3 θ 2π 0 dφ cos 2 φ |Y m l (θ, φ)| 2 = π dp p 4 |F nl (p)| 2 dθ sin 3 θ Y m l (θ, φ)| 2 (E64)