ORIGINAL ARTICLE: Bengyat, O., Di Biagio, A., Aspelmeyer, M., & Christodoulou, M. (2023). Gravity Mediated Entanglement between Oscillators as Quantum Superposition of Geometries. arXiv:2309.16312.
DOI: https://doi.org/10.48550/arXiv.2309.16312
Article Abstract: Protocols for observing gravity induced entanglement typically comprise the interaction of two particles prepared either in a superposition of two discrete paths, or in a continuously delocalized (harmonic oscillator) state of motion. An important open question has been whether these two different approaches allow to draw the same conclusions on the quantum nature of gravity. To answer this question, we analyse using the path - integral approach a setup that contains both features : a superposition of two highly delocalized center of mass states. We conclude that the two usual proto cols are of similar epistemological relevance. In both cases the appearance of entanglement, within linearized quantum gravity, is due to gravity being in a highly non - classical state: a superposition of distinct geometries.

Introduction

In this notebook we arrive at the result in Eq. (11), which is the generalized entanglement measure for the intermediate case of path superposition with non-negligible gaussian width. The calculation proceeds as follows: (1) defining the non-relativistic action, (2) finding the classical solutions of the system, (3) computing the gaussian integrals in Eq. (9), (4) calculating for this wavefunction the generalized entanglement measure [Swain et al, Phys. Rev. A 105, 052441 (2022)] defined in Eq. (28), to get Eq. (11).
<<VariationalMethods`​​​​$Assumptions={m>0,d>0,ω>0,ξ>0,t>0,ϕ>0,ϵ>0,Ω>0,ti>0,y1∈Reals,y2∈Reals,y1p∈Reals,y2p∈Reals,4ϵ^2ϕ<1,trp>=0,ξ1∈Reals,ξ2∈Reals,s1i∈Reals,s2i∈Reals,s1j∈Reals,s2j∈Reals,s1k∈Reals,s2k∈Reals,s1l∈Reals,s2l∈Reals};
Defining the dimensionless quantities ϵ, ϕ and Ω which fully define the problem .
ℏ=ϵ^2d^2mω;(*ϵ=α/disameasureofthe"quantumness"oftheparticles*)​​G=ϕℏdω/(m^2);(*ϕ=
2
Gm
/ℏωdisameasureoftheeffectofgravityonthesystem*)​​c=dω/(Ω);(*Ω=ωd/cisameasureoftheeffectofretardationinthesystem*)​​α=Sqrt[ℏ/(mω)];(*thewidthofthegaussian*)​​
A function calculating an integral of the form fPre Exp[fExp], where fPre is a polynomial and fExp a 2nd order polynomial. Taken from https://mathematica.stackexchange.com/a/6846.
gaussMoment[fPre_,fExp_,vars_]:=Module[{coeff,dist,ai,μ,norm},coeff=CoefficientArrays[fExp,vars,"Symmetric"True];​​ai=Inverse[2coeff[[3]]];​​μ=-ai.coeff[[2]];​​dist=MultinormalDistribution[μ,-ai];​​norm=1/PDF[dist,vars]/.Thread[varsμ];​​normExp[1/2coeff[[2]].μ+coeff[[1]]]Distribute@Expectation[fPre,varsdist]]
Functions replacing coefficients of given polynomials with single symbols. This makes Mathematica work much faster.
myMatrix[symb_,d1_,d2_,d3_,d4_]:=Table[Subscript[symb,i1,i2,i3,i4],{i1,d1},{i2,d2},{i3,d3},{i4,d4}];​​symbolize[poly_,vars2_,symb_]:=Module[{lst,flst,symbolized,dims,symbCoeffs},​​lst=CoefficientList[poly,vars2];​​dims=Dimensions[lst];​​symbCoeffs=myMatrix[symb,dims[[1]],dims[[2]],dims[[3]],dims[[4]]];​​Do[​​If[lst[[i1,i2,i3,i4]]==0,symbCoeffs[[i1,i2,i3,i4]]=0];​​,{i1,1,dims[[1]]},{i2,1,dims[[2]]},{i3,1,dims[[3]]},{i4,1,dims[[4]]}];​​symbolized=Internal`FromCoefficientList[symbCoeffs,vars2];​​symbolized​​];
A function wrapping gaussMoment, replacing coefficients with symbols before integrating.
myGauss[fPre_,fExp_,vars1_]:=Module[{a0exp,a0pre,aexp,apre,sympre,symexp,res},​​a0exp=CoefficientList[fExp,vars];​​symexp=symbolize[fExp,vars,aexp];​​a0pre=CoefficientList[fPre,vars];​​sympre=symbolize[fPre,vars,apre];​​res=gaussMoment[sympre,symexp,vars1]/.{aexp->a0exp,apre->a0pre,Subscript->Part};​​res​​]
Operations on complex numbers done in a custom way which for some reason works faster.
myRe[var_]:=Module[{var1,var1C,re,im},​​var1=var//ComplexExpand//Evaluate;​​var1C=var1//Conjugate//Refine;​​re=(var1+var1C)/2//Simplify;​​re​​]​​​​myIm[var_]:=Module[{var1,var1C,re,im},​​var1=var//ComplexExpand//Evaluate;​​var1C=var1//Conjugate//Refine;​​im=-(var1-var1C)/2//Simplify;​​im​​]​​​​myConj[var_]:=Module[{var1,var1C,re,im},​​var1=var//ComplexExpand//Evaluate;​​var1C=var1//Conjugate//Refine;​​var1C​​]
The series expansion to second order in ϵ and zeroth order in Ω. The second function is for quantities which will be taken with a square root later.
mySeries[expr_]:=Series[expr,{Ω,0,0},{ϵ,0,2}]//Normal;​​mySeriesSq[expr_]:=Series[expr,{Ω,0,0},{ϵ,0,4}]//Normal;

Action

We define action of the gravitational interaction: SGhb =
S
G
/ ℏ. This is the Newtonian potential, as we work in the non-relativistic limit in this calculation. The functions q1, q2 are the variables for the Lagrangian and action.
In[]:=
LG12=(Gm^2/2)/(d+αq1[ti]+αq2[ti])/ℏ//mySeries//Simplify;​​LGhb=LG12+(LG12/.{q1->q2,q2->q1})//mySeries//Simplify​​SGhb=Integrate[LGhb,{ti,0,t}];​​
Out[]=
ϕω(1-ϵ(q1[ti]+q2[ti])+
2
ϵ
2
(q1[ti]+q2[ti])
)
The free, uncoupled action: Sunc =
Σ
a=1,2
S
a
/ ℏ:
In[]:=
Lunc=(m(αq1'[ti])^2/2+m(αq2'[ti])^2/2)/ℏ//mySeries//Simplify​​Sunc=Integrate[Lunc,{ti,0,t}];
Out[]=
2
′
q1
[ti]
+
2
′
q2
[ti]
2ω
Now that we have defined the action, we proceed to finding the solution to the classical equations of motion, which we will use for the stationary-phase approximation in Eq. (9).

Classical Solutions

We find the classical trajectories of the particles for initial locations y1p, y2p =
y'
1
/α,
y'
2
/α
and final locations y1, y2 =
y
1
/α,
y
2
/α
. Here, the y’s are in units of α, the gaussian’s width. We do it by solving the system under the second order approximation in ϵ. The problem becomes two harmonic oscillators with some frequency trpω, which are coupled by another spring, which is gravity. The variable trp is just a bookkeeping variable which will be taken to zero to achieve the trajectories of the two otherwise free particles coupled by gravity.
In[]:=
L[ti]:=(LGhb+Lunc);​​L[ti]//Simplify
Out[]=
ϕω(1-ϵ(q1[ti]+q2[ti])+
2
ϵ
2
(q1[ti]+q2[ti])
)+
2
′
q1
[ti]
+
2
′
q2
[ti]
2ω
In[]:=
eqs=EulerEquations[L[ti],{q1[ti],q2[ti]},ti]//Expand(*theclassicalequationsofmotion*)
Out[]=
-ϵϕω+2
2
ϵ
ϕωq1[ti]+2
2
ϵ
ϕωq2[ti]-
′′
q1
[ti]
ω
0,-ϵϕω+2
2
ϵ
ϕωq1[ti]+2
2
ϵ
ϕωq2[ti]-
′′
q2
[ti]
ω
0
Their solution
sol=DSolve[{eqs[[1]],eqs[[2]],q1[0]==y1p,q2[0]==y2p,q1[t]==y1,q2[t]==y2},{q1[ti],q2[ti]},{ti}]//Simplify
x1[time_]:=q1[ti]/.sol[[1]]/.ti->time;​​x2[time_]:=q2[ti]/.sol[[1]]/.ti->time;​​x1[ti]​​x2[ti]
Out[]=
1
4-1+
4tϵ
ϕ
ω

tϵ
-2tiϵ
ϕ
ω

-
2tϵ
ϕ
ω

t(-1+2y1ϵ+2y2ϵ)+
2(t+2ti)ϵ
ϕ
ω

t(-1+2y1ϵ+2y2ϵ)+
4tϵ
ϕ
ω

t(-1+2y1pϵ+2y2pϵ)-
4tiϵ
ϕ
ω

t(-1+2y1pϵ+2y2pϵ)-
2tiϵ
ϕ
ω

(2ti(y1-y1p-y2+y2p)ϵ+t(1+2y1pϵ-2y2pϵ))+
2(2t+ti)ϵ
ϕ
ω

(2ti(y1-y1p-y2+y2p)ϵ+t(1+2y1pϵ-2y2pϵ))
Sanity check: initial and final locations
sanity check: check that the solution actually solves the classical equation

Wavefunction

Calculate the norm squared normsq of the state, integrating each separate exponent separately using our myGauss function:

Generalized Entanglement Measure

Each of the four multiplicands inside the integral above is slightly different in where is the ' and whether it is conjugate
Now sum over the 16 elements:
Finally, the generalized entanglement measure (Eq. 28) is:
Find a nicer representation (Eq. 11) with order 1 functions (Eq. 29):
Plot the functions:

CITE THIS NOTEBOOK

Gravity Mediated Entanglement between Oscillators as Quantum Superposition of Geometries​
by Ofek Bengyat​
Wolfram Community, STAFF PICKS, October 12, 2023
​https://community.wolfram.com/groups/-/m/t/3035123