Gravity-induced entanglement between two massive microscopic particles in curved spacetime: I. The Schwarzschild background

,


I. INTRODUCTION
One of the significant cornerstones in the history of physics is the establishment of quantum field theory, which has successfully described the interaction of all relativistic fields except gravity.Despite decades of substantial efforts, the quantization of the gravitational field remains a challenging endeavor, leaving the comprehensive theory of quantum gravity still elusive.A significant barrier lies in the absence of experimental evidence supporting the quantum aspects of gravity, despite several proposed experimental designs aimed at detecting quantum gravity phenomenology [1].
However, these previous studies assumed that massive particles were placed in an approximately flat spacetime (localized within the lab) and existed for short durations (seconds).In contrast, numerous astrophysical processes generate massive particles that interact over extended periods while traversing the universe, providing a natural setting for detecting the QGEM effect.This paper proposes an innovative scheme to demonstrate gravity-induced entanglement in more general curved spacetimes, involving smaller-scale particles (microscopic when compared to the mesoscopic particles utilized in the QGEM setting).Specifically, we investigate the generation of entanglement between two particles in a Schwarzschild background to universally and convincingly test the quantum gravity experiment.Assuming that the particle pairs move along geodesic paths in each instance, the separation between each pair of trajectories will change due to geodesic deviation [see Fig. 2].Consequently, the proper time will vary between each pair of trajectories due to alterations in spacetime geometry.Intuitively, the proper time of the closest particles, which have the shortest spacelike distance between them, will experience the most significant increase.Building upon [4], we compute the phase shift in each superposition state to determine the presence of entanglement.Remarkably, entanglement indeed emerges, illustrating the quantum gravity effect in general curved spacetime.We conduct a comprehensive analysis to explore the factors influencing the phase shift and their quantitative impact on entanglement.
However, a challenge remains in distinguishing whether the observed entanglement is generated by the quantum gravity effect of the particle pairs or by other processes during emission and propagation.To address this concern, we propose that QGEM during geodesic motion will exhibit a distinctive spectrum, as phase shifts occur across a series of geodesics.Analyzing the entangled patterns formed by various geodesics can assist in determining whether the entanglement originates from the gravitational field of nearby particles or from alternative sources.
The paper is structured as follows: In Sec.II, we provide a detailed description of our proposition for generating gravity-induced entanglement of microscopic massive particles in curved spacetime.In Sec.III, as an illustration, we consider a pair of particles moving within a globular clus-FIG.1.The QGEM experiment setup: In the initial phase of this experiment, two mesoscopic particles, each initially in a spin superposition state, are placed a short distance apart.Subsequently, the application of an inhomogeneous magnetic field prompts each particle to assume a spatially split state contingent on its spin configuration.This step ensures a spin-dependent spatial separation.Following this, a coherent superposition of states is maintained for a certain duration, while keeping the separation distance fixed within each branch.In the subsequent phase, the magnetic field is deactivated, and the spatially split states are realigned to regain their coherence.Finally, the experiment involves measuring spin correlations and calculating the entanglement witness.This analysis aims to determine if the system is indeed in an entangled state.Successful identification of entanglement would serve as confirmation of the quantum properties of gravity, consistent with the principles of entanglement theory.
ter with a Schwarzschild-like metric and explore the influence of initial conditions on gravity-induced entanglement.We analyze two galactic models with different mass profiles, namely the dual pseudo-isothermal elliptical density (dPIE) profiles and the Navarro-Frenk-White (NFW) profile.In Sec.IV, we examine a more realistic scenario where the particles' geodesic trajectories deviate from the galaxy's center.Additionally, to uncover more features of QGEM, we investigate the entanglement witness as a function of the particles' kinetic energy.The paper concludes in the final section.
Throughout this paper, we adopt the natural units system, c = G = 1, to simplify calculations.Physical quantities with units are provided in the SI system of units.

A. The scheme
In order to apply the QGEM scheme to astrophysics, its covariant description becomes crucial.The first covariant generalization of the QGEM scheme was undertaken by Christodoulou and Rovelli [4].They also highlighted that gravity-induced entanglement arises from the superposition of spacetime geometry, leading to distinct changes in proper time across the four branches.
To delve into specifics, the particle pair is initially prepared in a superposition state of spin and spatial position: where , and so forth.The notation |B⟩ represents the quantum state of the background gravitational field.
Subsequently, in accordance with the hypothetical quantum superposition of distinct spacetimes [4], the four branches will exhibit distinct time evolutions.The term |RL⟩, characterized by the shortest separation, will accrue the maximum phase due to the rapid growth of inherent time, as expressed by where m 0 signifies the particle's mass.Consequently, the central phase difference responsible for entanglement is represented by Upon recombining the two components of the superposition, the final state, accounting for an overall phase factor, can be articulated as Evidently, this state embodies entanglement of the spins of the two test masses, signifying that the gravitational field manifests as a quantum phenomenon.However, this covariant description [4] remains confined to a flat space-time background.Furthermore, within the QGEM setting, two masses are situated at the mesoscopic scale and are subjected to a Stern-Gerlach setting, leading to a superposition of two components where each particle occupies different positions.However, this operation becomes challenging for astrophysical sources located far from Earth.
To address these challenges, our approach deviates from the use of mesoscopic particles.In this study, we instead focus on microscopic massive particles, which had passed through a region with a magnetic field gradient during the propagation process in interstellar space and possessed a spatial superposition state.These microscopic particles freely fall within a generally curved spacetime, providing an opportunity to investigate gravity-induced entanglement.An overview of our general experimental setup is presented in Fig. 2.
Two identical microscopic particles, labeled A and B, both in the superposition state of two spatial positions as detailed in [2,3], traverse curved spacetime along their respective geodesic paths toward an Earth-based detector.Our main motivation is to employ a covariant methodology to comprehend the generation of entanglement throughout their journey.
In the reference frame of the moving particle, the particle's proper time is given by: τ = dτ .Let's now consider the scenario where another particle is in close proximity, moving alongside it but at a certain distance apart.The distance d (τ ) varies over time due to geodesic deviation.In this context, we assume that the gravitational attraction between the particles is significantly weaker than the tidal force, allowing d (τ ) to be primarily influenced by geodesic deviation.According to the equivalence principle in superposition spacetime [25] and the Newtonian limit approximation, the proper time for the two particles in each branch can be expressed as: where R ∼ λ c (λ c representing the Compton wavelength of the particle) denotes the radius of each body.It should be much smaller than the distance, R ≪ d(τ ), yet significantly larger than the Schwarzschild radius: R ≫ r s = 2m.The phase difference between each branch is entirely attributed to the term: resulting in a phase change that can now be expressed as: where m 0 denotes the static mass of the particles.Subsequently, we will primarily focus on calculating the phase change in one typical astrophysical scenario: a spherically symmetric background, such as a globular cluster.It is important to note that, in this experimental design, for this order of analysis, the moving particles should be massive.This is because a relatively stationary observer cannot perceive the gravitational effects of a massless particle, such as a photon.

B. Decoherence time
Our above scheme a prior assumes that the particles are in their spatial superposition states.However, in the complex background of the universe, various microscopic or mesoscopic particles, among which cosmic microwave background radiation (CMB) dominates, can potentially lead to decoherence.The decoherence rate of a neutral particle, initially in a superposition of two locations separated by ∆x and is immersed in a radiation heat bath at temperature T , can be described by the dipole approximation [26,27]: where ã represents the effective radius of the spherical particles.This formula holds when the radiation wavelength exceeds the particle radius and the superposition interval.For a superposition state with a particle mass of 10 −25 kg and an upper limit of the superposition interval of ∆x in the CMB background, the lower limit of the decoherence time is estimated to be 10 7 × m ∆x 2 s ∼ 10 19 × m ∆x 2 s.When this time scale greatly exceeds the Particle travel time, we are confidently assume that the superposition state is maintained throughout the entire particle movement.Of course, the scattering of impurity molecules in the universe and the thermal radiation of particles themselves will significantly increase the decoherence rate.

III. ENTANGLEMENT GENERATION BY GLOBULAR GALAXY
Our example focuses on entanglement induced by the gravitational field within a cluster of galaxies, teeming with neutral massive particles and possessing a size substantial enough to generate a pronounced entanglement effect.To provide a comprehensive exploration, we investigate two distinct galactic models with differing mass profiles.The first model employs the dual pseudo-isothermal elliptical density (dPIE) profiles [28,29], which have been demonstrated to aptly describe the mass distributions of the brightest cluster galaxies (BCGs) [30].These profiles find widespread use in lensing studies and deliver accurate fits to observed galaxies.The second model encompasses the Navarro-Frenk-White (NFW) profile [31], derived from N-body simulations, and widely employed for simulating dark matter (DM) halos within the ΛCDM universe.

A. dPIE model
In this subsection, let us envision an isolated galaxy cluster within the universe.For the sake of simplicity, we will exclusively consider the brightest cluster galaxy (BCG) component of the cluster.As previously indicated, the spherical dual pseudo-isothermal elliptical (dPIE) profiles prove especially fitting for characterizing the mass distributions of BCGs.These profiles are defined by their 3D-density [28][29][30]: where r is the distance from the center of the mass distribution, r core is the core radii and and r cut is the truncation radii with r cut > r core .While ρ 0 is the central density, which is related to the 1D-central velocity dispersion, σ 0 , by [28] ρ 0 = σ 0 2πG Integrating Eq. ( 9), one can obtain the total mass m(r) enclosed by a sphere of radius r when r < R, where R is the total radius of the BCG.Hence, the total mass of the galaxy is when r > R. The complete set of parameters utilized in the model adopted for this paper is enumerated in Table I.
The parameters of the BCG model are grounded on the fitting conducted on the A383 cluster, as demonstrated in [30].
The parameter values for NFW, which will be employed in the subsequent subsection, also derive from the A383 dataset.However, to enable a meaningful comparison between the two models, suitable adjustments have been applied to ensure that both models possess identical radius and mass values.With this mass profile in consideration, we are now prepared to present the metric for the background.Particularly, beyond the BCG, i.e., for r > R, the metric can be succinctly described by the Schwarzschild metric, with the mass specified in Equation ( 12) [32].
(13) While for the interior of the BCG (i.e., r < R), the metric is given by [see Appendix A] Here A(r) is obtained by integrating x 2 dx and the explicit expression is given in Eq. (A2), and B(r) is given by For convenience in the following calculations, let us unify these two cases into one form like Eq. ( 14), i.e., From now on, let us focus on one specific trajectory as shown in Fig. 3.One of the particles starts at r 0 , moving toward the center of the cluster and arrive the observer at r f finally.The whole trajectory follows a geodesic line through the center of the cluster.
It turns out the tetrad formalism is particularly convenient to discuss the geodesics.As to the present case, we construct the following four unity vectors, so that they form an orthogonal base where Recalling the geodesic equations in terms of four-velocity u µ (τ ) is For radial geodesics passing through the center (the origin) as shown in Fig. 3, the geodesic equations, after substituting Eqs. ( 14) and ( 17) into Eq.( 20), reduce to where the prime denotes differentiation with respect to r, while the dot represents differentiation with respect to τ .Note that the full trajectory can be split into four stages as shown in Fig. 3.The above equations (21) are only used for the first two stages.However, using the spherically symmetric nature of the spacetime, the third and forth stages can be viewed as a reverse of the stages 2 and 1, respectively.Therefore, equations ( 21) can be used by redefining Our first goal is to solve these equations (21), which are difficult to solve analytically.We therefore find the numerical solutions.Without loss of generality, for later numerical computations, we assume that the particle has static mass m 0 = 10 −25 kg, which is a typical mass of a microscopic particle (with the Compton wavelength λ c ∼ 10 −15 m).
Fig. 4 shows a set of particle trajectories.The initial value for radial coordinate r 0 in our numerical calculations is taken to be 4.5 × 10 16 (in natural unit 1 ), while the initial values for ṙ(τ = 0) and ṫ(τ = 0) can be transformed to the initial values 1 Translated to the International System of Units (SI system), r 0 = 1.35 × 10 22 km ∼ 4.375 × 10 8 pc.
of four velocity as where and v 0 is the initial velocity of the particle.

Fig. (4a)
illustrates that the value of r decreases almost linearly with τ , implying the particle's propagation from the source to the cluster.The steeper slope of the curve corresponds to a larger initial velocity.This correlation is reasonable since particles with higher v 0 will outpace those with lower v 0 .In Fig. (4b), the curves of t vs. τ largely overlap across the four different v 0 values, suggesting insensitivity to the initial velocity's magnitude.
Taking into account that the spacelike interval between two adjacent branches must exceed the Compton wavelength, we provide distinct initial conditions for the following scenarios.Now, let's proceed to calculate the geodesic deviation vectors, from which values of d(τ ) can be deduced-this term plays a pivotal role in evaluating δϕ as shown in Equation (7).These vectors are determined through the solution of geodesic deviation equations in the tetrad formalism.Unlike the fixed tetrad used in Equation ( 17) for geodesics, the tetrad here must be parallelly transported along the geodesics [33].In other words, ẽ(a) µ ;ν v ν = 0, where v ν is the tangent vector of the geodesics.Consequently, the axes' orientations remain fixed and non-rotating, as ascertained by local dynamical experiments.It turns out that the following selection of tetrads is appropriate ẽ( 1) Note that in order for the above tetrad satisfying the orthogonal normalized condition (19), one extra condition should be imposed The timelike basis ẽ(0) µ is the geodesic four-velocity vector and the remaining three are spacelike basis pointing in different direction.All the four basis compose an orthogonally translational basis.In this tetrad basis the metric becomes just Minkowski metric.
Then according to the geodesic deviation equation in the above frames, we have where w (a) is geodesic deviation vector in tetrad form and Substituting the metric ( 14) and the tetrad ( 24) and ( 27) into above definition, one finds that the nonvanishing components of k (a)  (b) are Without loss of generality, we assume that w µ is orthogonal to ẽ(0) µ , such that we have w (0) = 0. Due to the spherical symmetry of the background, the second and third components of w (a) are simply symmetric angular components.Therefore, we can also choose suitable frame such that the initial value of w (3) is vanishing.Note that Eq. ( 28) is ho-mogenous, hence if ŵ(a) is a solution of this equation, so is κ ŵ(a) , where κ is an arbitrary nonzero constant.Hence, once we initially have vanishing w (3) and ẇ(3) , it will always be.
In the end we are left with two nonzero components w (1) and w (2) .
The total spacelike distance d(τ ) between the pair of particles is where ˜means the distance in the stages 3 and 4. ∆l is a constant and can be absorbed into w (a) in our numerical computations.Integrating over all four stages, one leads to the change of the proper time Once δτ is obtained, the change in phase can be known through As a consequence, the key step is to get the values of w (i) and w(i) (i = 1, 2), which can be obtained by solving the geodesic deviation equations ( 28), (31) and (32).
In what follows let us turn to numerically solve the geodesic deviation equations ( 28), (31) and (32) under the metric ansatz Eqs. ( 14)-( 16) and (A2), and obtain δϕ for appropriate initial conditions.Totally the initial conditions include: the initial geodesic deviations w (1) 0 and w (2) 0 , the initial geodesic deviation velocities ẇ(1) 0 and ẇ(2) 0 measured in tetrad basis (24), the initial radial coordinates r 0 , and lastly the initial values for ṙ(τ = 0) and ṫ(τ = 0), which, again, can be transformed to the initial values of four velocity via (23) by setting suitable initial velocity v 0 .Totally there are six initial conditions to be assigned.However, for simplicity, in our following numerical simulations we will assume that w The full set of initial values which are adopted in the numerical simulations can be found in Table II.Note that throughout the paper, z represents the redshift of the source 2 .It is not difficult to generalize the simulations to the cases where w 2 The redshift z can be calculated by using where Ω Λ , Ω M and Ω R are the density parameters for dark energy, matter and radiation, respectively.H 0 is the Hubble constant.According  Fig. 5 displays geodesic deviation for various initial values, revealing significant changes in all space-like separations d(τ ) near the center of the galaxy.Notably, the magnitude of d directly affects the extent of change when approaching the galaxy center.In some cases, like curve S1 in Fig. 5a, where d(τ ) is very small near the center, the turning point might not be clearly represented on the graph (indicated by the black dashed line in Fig. 5 as the particle reaches the center of the cluster).For particle pairs with a negative initial geodesic deviation velocity (particles initially moving towards each other), they initially converge and subsequently diverge.This behavior is consistently captured in our numerito Plank 2018 data [34]: Ω Λ = 0.685, Ωm = 0.315, Ω R = 0 and H 0 = 67.4km/s/Mpc.Without loss of the generality, in the paper we take r f = 4.5 × 10 16 for simplicity.cal results, as evidenced by the curve of set S3 in Fig. 5, which first declines and nearly touches the zero point before passing through the center of the cluster, and then rises afterward.Phase changes δϕ in different initial conditions are shown in Fig. 6.In each subgraph, the shade of color represents the magnitude of the phase change varies with two of four initial values, with the remaining two values keep fixed.
From the left six pictures, it is obvious that the phase change is sensitive to the initial geodesic deviation velocity ẇ0 no matter whether it is positive or negative.More specially, from figures in the leftmost panel one can see that the phase change increases at the beginning with negative deviation velocity − ẇ0 and then decreases.This is contrary to that of d(τ ) as shown in Fig. 5a.This is reasonable because the separation distance d(τ ) appearing in the denominator of Eq. (34).While figures in the second column show that δϕ decreases with initial deviation velocity ẇ0 monotonously, which is also consistent with the results as shown in Fig. 5.In both cases, δϕ increases monotonously with the redshift z, although it is relatively mild for the negative initial deviation velocity case.In contrast, the phase change decreases monotonously both with the initial particle velocity v 0 and initial deviation distance w 0 no matter whether the initial geodesic deviation velocity is positive or negative.
On the other hand, three figures in the bottom right corner corresponds to the cases where the initial deviation velocity is 0, leaving dependence of δϕ on z, v 0 and w 0 , respectively.Specifically, the third column shows that for fixed z, δϕ decreases with both v 0 and w 0 monotonously.In contrast, it increases with z especially obvious for fixed v 0 .This is because as z increases, the proper time of the whole geodesic motion also increases, but the geodesic deviation is basically unaffected.The rightmost figure shows that for fixed z and v 0 , the phase change is inversely proportional to the initial geodesic deviation w 0 .This is because the geodesic deviation equations are a set of homogeneous second-order differential equations and the spacelike interval at every moment is proportional to its initial value w 0 .

B. NFW model
In this subsection let us consider the NFW model, which can be useful to observe the effects of the DM halos in the ΛCDM universe.The density profile and mass of the NFW model is [30] ρ NFW (r) = ρ 0 where ρ 0 is the scaling density and r s is the scaling radius.These two parameters have relationship with the Virial mass of the halo through M 200 ∝ ρ 0 r 3 s and are related to the concentration parameter through c = r 200 /r s .In order to compare with the dPIE model, in this paper we set ρ 0 such that the total mass of the cluster is the same as the one of the dPIE model.The full list of the parameters can be found in Table I.Following the same procedures, one finds the metric parameters A(r), B(r) in ( 14) is FIG. 6. Contour plot of phase changes δϕ for initial value set S1 in Table II.In these plots, we let two of four initial conditions vary in each subgraph, while the other two are fixed.Different color denotes different values of the phase change.The following notations w Following what we did in the last subsection, one can obtain phase change δϕ for different initial conditions.Fig. 7 depicts the phase change δϕ as a function of initial redshift, considering other initial values given in Table II.Notably, the plot exhibits a nearly linear growth of the phase change with increasing initial z.This behavior is attributed to the increase in the proper time of motion before passing through the center of the galaxy cluster, which serves as the integral variable, while the geodesic deviation after passing through the center remains almost unchanged.Consequently, the curve closely resembles a straight line.
In Fig. 8, we plot the phase change δϕ as a function of initial radial velocity v 0 .One can see from this plot that the phase change decreases monotonously as increasing v 0 .In addition, for small v 0 , it decreases fast as v 0 increases and gradually the decrease becomes mildly.Finally it approaches zero gradually for large v 0 .This is can be explained in the following way: the particle's total proper time integrated over the trajectory decreases precipitously as initial radial velocity v 0 goes from zero to non-zero.This has been already verified in the dPIE model as shown in Fig. 5b.
As we can see from Fig. 9, the phase change is just inversely proportional to w 0 as expected.Fig. 10 is a plot of δϕ as a function of initial geodesic deviation velocity ẇ0 of the particles 3 .The initial values of other parameters are given in Table II.From Fig. 10a, we find that for smaller ẇ0 , δϕ decreases slowly with increasing ẇ0 , then suddenly for some point ( ẇ0 ≈ 10 −37 ), it decrease quickly and then tends to be nearly flat.This behavior can be roughly explained in the following way: Eq. ( 28) can be converted to d 2 w (a) dτ 2 = −k (a) (b) w (b) , the r.h.s of the equation looks like a "geodesic deviation acceleration" which is proportional to k (a) (b) and w (b) .Meanwhile, comparing with other regions, k (a)  (b) has dominant effects near the center of the cluster.As a consequence, the particle pair with a larger w (b) will acquire an overwhelming geodesic deviation acceleration as they pass 3 Again, for simplicity, we have assumed that ẇ(1)  II.
the center of the cluster.The separation distance d(τ ) between them will continue to increase after this process.So when the absolute value of the geodesic deviation velocity reaches a certain value, w (b) (or d(τ )) will be large enough and the phase change will drop sharply as it is inverse proportional to d(τ ) as shown in Eqs.(34) and (35).The negative initial geodesic deviation velocity case, how-  II.ever, is very different as shown in Fig. 10b.δϕ of this case increases slowly with increasing − ẇ0 at the beginning, and then suddenly increases to a maximum at some point ( ẇ0 ≈ 10 −37 ), followed by a quick decrease to lower value.Finally it descends gradually to some finite value.Recalling a similar behavior existing in S3 of Fig. 5a, where d(τ ) firstly goes down quickly as increasing − ẇ0 and then grows up vastly.We know that this is because for the particle pairs with a negative initial geodesic deviation velocity, at the beginning the particle pairs approach each other and then they move away from each other.Similarly, in the present case, the particles with a negative initial geodesic deviation velocity approach each other firstly and later will move away from each other, so that d(τ ) will decrease quickly at the beginning and then increase fast, leading to inverse behavior of δϕ as d(τ ) appearing in the denominator of Eq. ( 34).
In the previous simulations, we let the total mass M of the cluster fix.In order to see the influence of M on the entanglement phase, it is useful to calculate the phase change with different mass in dPIE and NFW model.We first plot phase change as a function of M in the dPIE model as shown in Fig. 11a.Then we define a ∆ϕ = δϕ dP IE − δϕ N F W , where δϕ dP IE and δϕ N F W denote, respectively, the phase change for dPIE and NFW models.We then plot ∆ϕ as a function of M as shown in Fig. 11b.We can find from Fig. 11a that phase change decreases faster and faster as M increases.This is because the stronger the gravitational field, the shorter the proper time, the phase change decreases with galaxy's mass.Then we can see from Fig. 11b that the specific model also act on the phase change in a certain way.In the smaller mass range NFW model is more conducive to the formation of phase change while in the bigger mass range dPIE is better at inducing phase change.
In summary, compared with dPIE model, the relationship between the phase change and each initial value exhibits the similar characteristics.Generally speaking, larger z, smaller initial radial velocity v 0 , and appropriate negative initial deviation velocity lead to bigger δϕ and (possibly) more significant effects of entanglement.

IV. OFFSET GEODESICS AND CHARACTERISTIC SPECTRUM
In the preceding discussion, we made the assumption that all particles follow geodesics passing through the cluster's center.However, in reality, geodesics may deviate from the cluster's center in various ways.Moreover, from an experimental standpoint, distinguishing whether observed entanglement arises from the quantum gravity effect of particle pairs or from their emission at the source presents a significant challenge.Therefore, it becomes crucial to discern whether gravity induces entanglement or if other physical processes are at play.
To tackle this challenge, we propose exploring the characteristic spectrum of QGEM during geodesic motion, which manifests as phase changes along a series of geodesics.Through analysis of the entangled patterns formed by dif-  I.The initial values of the other parameters are given by S1 in Table II.
ferent geodesics, we can infer the entanglement's origin-whether it arises from the gravitational field of nearby particles or from alternative sources.
Generally speaking, there are two possible geodesics, contingent upon our knowledge of the particles' source location, as illustrated in Fig. 12.In what follows, we would like to discuss possible characteristic spectrum of the QGEM in these two cases, respectively.

A. The particles' source location is given
In the first geodesic case, we consider a globular cluster with a density profile described by the dPIE model [30].The parameter values of the model are chosen the same as in Table I, and the particle's mass is assumed, again, to be 10 −25 kg.For the sake of simplicity, in this subsection, we only consider the set S1 from Table II as the initial parameter values.
At the source point, we use the same orthogonal frame as in (17) and (18).The four-velocity is expressed as: where γ = and v 0 is the particle's initial emission velocity.The four-velocity mentioned above is connected to the deflection angle ξ, with ξ being assigned to values in the range of ξ ∈ [0, π) for the initial four-velocity.For particles that can reach the Earth, ξ and v 0 are not independent.Namely, for each ξ, there is a unique value of v 0 .To simplify our analysis, we will focus on the scenario where the entire trajectories lie in the same plane.Under these initial conditions, the geodesics will deviate from the center of the cluster.Notably, when ξ = 0, the geodesics that pass through the galaxy's center will be recovered.We also build a set of parallelly transported tetrads and assume the initial geodesic deviation vector is: where ẽ(∥) is the basis parallel to the connecting line between the Earth and the galaxy cluster , ẽ(⊥) is the basis perpendicular to ẽ(∥) and four-velocity.The initial geodesic deviation component w ∥ = 1 3 × 10 −16 and w ⊥ = 0. Now let us delve into the QGEM effect in this scenario.As usual, we can denote positive helicity as ↑ and negative helicity as ↓.Therefore, the adjacent propagating particles are in a superposition state of 1 change induced by gravity is again described by (7).Following the same procedures given in section III, one can calculate the phase change under different initial conditions.Compared to the initial conditions, particle's energy is an important observable.In observer's frame, each particle emitted at a different angle will have a specific speed, and the Earth's observer will measure a corresponding kinetic energy.In this static spherically symmetric space-time, there is a conserved quantity along geodesics, denoted as C [4], is given by 4 : where ξ µ represents the timelike killing vector (∂t) µ , u µ is the particle's four-velocity (38), and E ob is the energy measured by a local stationary observer at infinity (e.g., an observer on Earth).We use the kinetic energy of the particles measured by observer on Earth to label each particle.The observed kinetic energy for particles emitted from different initial angles is calculated as follows: where, r earth represents the coordinate distance between the center of the Earth and the center of the galaxy cluster, and r 0 represents the coordinate distance between location of particle emission and the center of the galaxy cluster.
The result is shown in Fig. 13.It shows that as the emission angle ξ increases, the kinetic energy of particles reaching Earth decreases, and the intrinsic length of the geodesic line increases, leading to a progressively larger entanglement phase.The irregular jumps on the curve are caused by the instability of numerical simulation.
Phase change is not directly observable.On the other hand, the entanglement witness W is often utilized as an experimental indicator to detect entanglement formation.The definition of entanglement witness is 4 This is because: for a given Killing vector ξ µ and a tangent vector u µ = (∂τ ) µ of a geodesic γ(τ ), the identity holds: u µ ∇µ(ξν u ν ) = 0.When it's greater than 1, we can infer that there is entanglement between the two particles.The variation of W with respect to E kin is plotted in Fig. 14.From this figure, we observe that the entanglement witness oscillates rapidly with energy, but it exhibits different oscillating frequencies in various energy segments.Specifically, it oscillates relatively fast in the region of low energy, which coincides with the fact that the phase increases faster with an increase in the emission angle.The vertical axis represents the entanglement phase δϕ formed during the propagation process, while the horizontal axis corresponds to the kinetic energy E kin measured by the Earth observer after the particles with different emission parameters reach the Earth.

B. The particles' source location is unknown
In this case, since we lack information about the location of the particles' source, we cannot use the deflection angle to span the geodesics.Instead, we utilize h, which represents the initial vertical distance from horizontal geodesics, to effectively span the geodesics, as illustrated in Fig. 12b.Each particle exhibits a distinct starting h value and, when coupled with an appropriate initial four-velocity, can be successfully detected by a probe on Earth.The four-velocity can be expressed as follows: where e ∥ µ is the unit space-like base formed by the linear combination of e (1) µ and e (2) µ , making it parallel to the connecting line between the Earth and the center of the galaxy cluster.Namely, In this case, the kinetic energy observed by the observer on Earth is determined by and each geodesic to the Earth corresponds to a specific value of v 0 and γ.The characteristic entangled spectral lines are shown in Figure 15 and Fig. 16.Similarly, the entanglement phase drops rapidly with increasing energy, and the entanglement witness oscillates slowly at high energy and rapidly at low energy.This behavior indicates that at greater emission heights, the phase increases faster.The plot of entanglement witness as a function of the emission height h can show more details.From Fig. 17, we see that for dPIE models with different center densities, the segment of the entanglement curve where FIG. 17. Variation of entanglement witness with particle emission height for dPIE models with different center density parameters.
Here, we assume that the value range of h is 0.05R to 0.2R.Other galaxy parameters and kinematic parameters remain unchanged.
the entanglement exceeds 1 varies.Additionally, the entanglement witness with greater center density oscillates more slowly.In other words, a higher center density ρ 0 of a galaxy will impede the generation of entanglement between particles.A possible explanation is that, under the condition of the same emission height, particles take less proper time to reach Earth when passing through more massive galaxies, which hinders entanglement formation.In summary, in both cases, the entanglement phase varies monotonously with the initial geodesic emission parameters, ξ and h, causing the entanglement witness W to oscillate with the initial parameters (or the kinetic energy), as shown in Figs.14, 16 and 17.This characteristic could be a significant index to distinguish the origin of the observed entanglement.More specifically, if the particles' entanglement is native (formed during particle emission), then the entanglement witness will be randomly distributed with energy and will not exhibit the quasi-periodic oscillation behavior mentioned above.In other words, only entanglement induced by quantum gravity effects can result in W exhibiting this quasi-periodic behavior.Thus, the identification of these entanglement features in experimental observations can provide confidence in concluding that the entanglement is induced by quantum gravity along the geodesic, rather than being formed during their emission.
The above conclusions are also valid for the NFW model.Fig. 18 and Fig. 19 show the entanglement witness as a function of kinetic energy, using the initial NFW parameters given in Tables I and II.These plots explicitly exhibit similar quasi- oscillation behavior as observed in the dPIE model in BCG.

V. CONCLUSIONS AND DISCUSSIONS
In the realm of quantum gravity, the challenge lies not in the absence of a complete mathematical physical theory, but rather in the scarcity of experimental approaches to connect theory and reality.In light of recent advancements in quantum technologies, endeavors have been made to elucidate the quantum nature of gravity.Among these endeavors, the quantum gravity induced entanglement of masses (QGEM) proposal has garnered significant attention [2][3][4]35].In this paper, we extend the QGEM experiments to include curved spacetimes.More specifically, we consider the QGEM for a pair of particles traveling along their geodesics in a galaxy with Schwarzschild metric as the spacetime background.We find that particle pairs readily become entangled at larger radial coordinates with appropriate small initial radial velocities.
By investigating the relations between entanglement witness and the kinetic energy observed by the observer (determined by the initial values of the model parameters), we find that there has a characteristic spectrum of QGEM.This provides us a way to distinguish whether the observed entanglement arises from the quantum gravity effect of particle pairs or from other process (e.g, the emission stage at the source as suggested by the Hawking radiation of black holes).
In addition, execution of the proposals outlined in this paper will demonstrate a more comprehensive quantum gravity effect in extensive spacetime, surpassing the limitations of local experiments carried out in labs [2,3].Notably, our scheme accommodates much lighter particle masses and significantly larger particle spacings.These particles may traverse hundreds of millions of years of geodesic movement before detection, an impractical feat on Earth.Additionally, direct detection of particles from the universe simplifies the experiment's preparation process.Microscopic particles may spontaneously exist in a superposition state somewhere in the universe, coinciding with our computational hypothesis mentioned above.
Furthermore, our astronomical scheme, based on the hypothesis of the equivalence principle in superposition spacetime [25], hopefully serves as an alternative validation for this extended equivalence principle.
Admittedly, there are some technical details need improve in the future.The geodesic deviation equation serves as a firstorder approximation formula for describing spacelike distance variations between particles.We solely consider the background spacetime metric while neglecting the backreactions of the nearby particles.Numerical calculations introduce method error and truncation error in the numerical simulation.Additionally, our idealized assumptions may not fully capture the complexities of the real universe.
In future research, we aim to employ more sophisticated calculation methods to account for the gravitational force of nearby particles.Additionally, we will explore more complex geodesic trajectories beyond radial movement to detect as many entangled particles as possible.Furthermore, in cosmic spacetime, two particles may be separated by considerable spacelike distances.As highlighted in [35], we will seek new methods to extract spacelike entanglement and investigate how gravity induces spacelike entanglement in curved spacetime.

FIG. 2 .
FIG. 2. Illustration of particle trajectories.The particles labeled as A and B are identical and microscopic in nature.They exist in a superposition state across two distinct spatial positions.The solid line delineates the geodesic trajectory of the particles within the curved spacetime.

FIG. 3 .
FIG. 3. Sketch of particle trajectory in galaxy cluster.It includes four different stages.The first stage refers to the trajectory from the starting point (r = r0, φ = φ0) to the boundary of the cluster (r = R, φ = φ0).The second stage involves the boundary r = R to the center of the cluster (r = 0).The third stage is the reverse of the second one, namely, from (r = 0) to (r = R, φ = φ0 + π).The last stage, which can be regarded as the reverse of the first stage, involves the geodesics from (r = R, φ = φ0 + π) to (r = r f , φ = φ0 + π), the observer.

FIG. 4 .
FIG. 4. Particle trajectories for different initial particle velocity v0with fixed r0 = r f ∼ 4.5 × 10 16 .(a) Relations between radial distance r and the proper time τ .Value of r decreasing with τ indicates that the particle is propagating from the source to the cluster.Negative values of r represent that the particle has passed through the center of the cluster.(b) Relations between temporal coordinate t and the proper time τ .They are insensitive to the initial velocity.

1
Note that this is in natural unit system.In SI system, v0, w by c = 3 × 10 8 m/s.For example, for S1, we have v0 = 3 × 10 4 m/s, w

FIG. 5 .
FIG.5.Spacelike separation d(τ ) between particle pairs as a function of initial conditions.The black dashed line denotes the moment when the particle reaches the center of the cluster.(a) Overview of particles start for different initial conditions as listed in TableII.(b) Full view of initial value set S5. Particle's proper time integrated over the trajectory of this case is much longer than those of others.

FIG. 7 .FIG. 8 .
FIG. 7. The phase change δϕ as a function of initial radial position (redshift) of the particles.The initial values of the other parameters are given in TableII.

FIG. 9 .
FIG.9.The phase change δϕ as a function of initial geodesic deviation w0 of the particles.The initial values of the other parameters are given in TableII.

FIG. 10 .FIG. 11 .
FIG. 10.The phase change δϕ as a function of initial geodesic deviation velocity ẇ0 of the particles.The initial values of the other parameters are given in TableII.

FIG. 12 .
FIG. 12. Two series geodesics offset from the center of the galaxy.(a) The particles' source location is given.(b) The particles' source location is unknown.

FIG. 13 .FIG. 14 .
FIG.13.Plot of entangled phase and observed particle energy.The vertical axis represents the entanglement phase δϕ, which is formed in the propagation process.On the horizontal axis, we have the kinetic energy Ekin measured by Earth observers after the particles with different emission parameters have reached the Earth.

FIG. 15 .
FIG.15.Entangled phase as a function of observed particle energy.The vertical axis represents the entanglement phase δϕ formed during the propagation process, while the horizontal axis corresponds to the kinetic energy E kin measured by the Earth observer after the particles with different emission parameters reach the Earth.

FIG. 16 .
FIG.16.Entanglement witness as a function of kinetic energy.

FIG. 18 .FIG. 19 .
FIG.18.Entanglement witness as a function of kinetic energy in the NFW model as the location of the particles' source is given.

TABLE I .
Parameter choices of two different galaxy models

TABLE II
. Sets of initial values for numerical simulationsSets redshift z velocity v0 w