PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
Nume ical simula ion o la ge-scale nonlinea open quan um mechanics
M. Roda-Llo des ,1,2D. Candoli ,1,2P. T. G ochowski ,1,2,3A. Rie a-Campeny ,1,2T. Ag enius ,1,2J. J. Ga cía-Ripoll ,4
C. Gonzalez-Balles e o ,1,2and O. Rome o-Isa 1,2
1Ins i u e o Quan um Op ics and Quan um In o ma ion o he Aus ian Academy o Sciences, 6020 Innsb uck, Aus ia
2Ins i u e o Theo e ical Physics, Uni e si y o Innsb uck, 6020 Innsb uck, Aus ia
3Cen e o Theo e ical Physics, Polish Academy o Sciences, Aleja Lo ników 32/46, 02-668 Wa saw, Poland
4Ins i u o de Física Fundamen al IFF-CSIC, Calle Se ano 113b, Mad id 28006, Spain
(Recei ed 19 June 2023; accep ed 8 Feb ua y 2024; published 8 Ma ch 2024)
We in oduce a me hod o sol e nonlinea open quan um dynamics o a pa icle in si ua ions whe e i s s a e
unde goes signi ican expansion in phase space while gene a ing small quan um ea u es a he phase-space
Planck scale. Ou app oach in ol es simula ing wo s eps. Fi s , we ans o m he Wigne unc ion in o a ime-
dependen ame ha le e ages in o ma ion om he classical ajec o y o e icien ly ep esen he quan um
s a e in phase space. Nex , we simula e he dynamics in his ame using a nume ical me hod ha implemen s
his ime-dependen nonlinea change o a iables. To demons a e he capabili ies o ou me hod, we examine
he open quan um dynamics o a pa icle e ol ing in a one-dimensional weak qua ic po en ial a e ini ially
being g ound-s a e cooled in a igh ha monic po en ial. This app oach is pa icula ly ele an o ongoing e o s
o design, op imize, and unde s and expe imen s a ge ing he p epa a ion o mac oscopic quan um supe posi ion
s a es o massi e pa icles h ough nonlinea quan um dynamics.
DOI: 10.1103/PhysRe Resea ch.6.013262
I. INTRODUCTION
The ield o le i odynamics [1], which ocuses on le i a ion
and con ol o mic oobjec s in acuum, allows us o s udy he
cen e -o -mass mo ional dynamics o a pa icle in a highly
isola ed en i onmen . Since he mechanical po en ial in which
he pa icle mo es can be con olled bo h dynamically [2–4]
and s a ically [5,6], le i a ed pa icles o e a unique pla -
o m o s udy nonlinea conse a i e mechanics. Fu he mo e,
he cen e -o -mass he mal ene gy can be emo ed, ei he
ia ac i e o passi e eedback, o he ul ima e limi whe e
only quan um luc ua ions a e p esen [7–13]. Cen e -o -mass
g ound-s a e cooling and he con ol o he mechanical po en-
ial open he possibili y o s udy nonlinea quan um mechanics
wi h a mic osolid con aining billions o a oms [6,14]. To
design, op imize, and unde s and expe imen ally easible p o-
ocols in ol ing nonlinea quan um mechanics, i is c ucial
o ha e a eliable nume ical ool ha allows us o e icien ly
simula e he dynamics while accoun ing o sou ces o noise
and decohe ence. In his pape , we p o ide such a ool in he
pa icula ly ele an and challenging scena io o mul iscale
dynamics induced by cen e -o -mass cooled massi e pa icles
e ol ing in wide nonha monic po en ials.
Mo e speci ically, he cen e -o -mass mo ion o cooled
mic opa icles exhibi s minu e luc ua ions (i.e., ze o-
poin mo ion), smalle han he size o a single a om.
Published by he Ame ican Physical Socie y unde he e ms o he
C ea i e Commons A ibu ion 4.0 In e na ional license. Fu he
dis ibu ion o his wo k mus main ain a ibu ion o he au ho (s)
and he published a icle’s i le, jou nal ci a ion, and DOI.
Expe imen ally easible nonha monic po en ials a e wide
han ze o-poin mo ion leng h scales, ha is, he dis ance
be ween classical u ning poin s is o de s o magni ude la ge
han he ze o-poin leng h scale. Hence, he dynamics ig-
ge ed in hose nonha monic po en ial will gene a e la ge
phase-space expansions. This expansi e dynamics will e en-
ually ac i a e he nonha monici ies in he po en ial, such as
a u ning poin s, which in he case o cohe en dynamics
can c ea e phase-space s uc u es a o e en below he Planck
scale [15]. This mul iscale phase-space dynamics o he pa -
icle’s cen e -o -mass s a e will be s udied h ough he ime
e olu ion o he co esponding Wigne unc ion. The use o
he Wigne unc ion is ad an ageous as i enables o inco po-
a e sou ces o noise and decohe ence (i.e., open dynamics)
while also clea ly iden i ying quan um ea u es [16] (e.g.,
h ough nega i e alues in he Wigne unc ion). The amoun
o compu a ional esou ces needed by nume ical me hods
based on ixed g ids [17] scales quickly wi h he amoun
o expansion in he dynamics. Thus, desc ibing he scena io
o in e es , wi h expansions o se e al o de s o magni ude,
becomes un easible using such me hods. The e o e, a mo e
e icien ep esen a ion o he s a e o his speci ic dynam-
ics is necessa y. We p opose using a dynamical g id me hod
whe e we desc ibe he s a e in a ime-dependen phase-space
ame ela ed o he o iginal ame acco ding o he classical
ajec o ies dic a ed by he nonha monic po en ial. As we
show below, nume ically simula ing he dynamics in such
a ame using a ixed g id is equi alen o a physically in-
o med adap i e g id ha places he g id poin s whe e hey a e
mos ele an , he eby imp o ing compu a ional e iciency.
This app oach has p o en in aluable in he design, op imiza-
ion, and unde s anding o a ecen p oposal o gene a ing
2643-1564/2024/6(1)/013262(10) 013262-1 Published by he Ame ican Physical Socie y
M. RODA-LLORDES e al. PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
mac oscopic quan um supe posi ions o a nanopa icle
h ough he nonlinea quan um mechanics induced in a wide
double-well po en ial [6].
This pape is s uc u ed as ollows: In Sec. II, we p esen
he heo e ical amewo k o ou me hod, including he ime-
dependen change o a iables leading o he ime-dependen
phase-space g id. In Sec. III, and in a dedica ed Appendix, we
de ail ou nume ical implemen a ion using ini e di e ences
and classical ajec o y p opaga ion. We hen examine he
dynamics in weak qua ic po en ials as an example o la ge
expansions wi h Planck-scale quan um ea u es in Sec. IV.
Finally, we conclude wi h ou inal ema ks and ou look in
Sec. V.
II. WIGNER FUNCTION DYNAMICS
IN THE LIOUVILLE FRAME
We conside a pa icle wi h mass me ol ing in a one-
dimensional po en ial U(x) in he p esence o noise. We
desc ibe he s a e o he pa icle h ough i s Wigne unc ion
W(x,p, ). The equa ion o mo ion o he Wigne unc ion is
gi en by
∂W(x,p, )
∂ =(Lc+Lq+Ln)W(x,p, ).(1)
The i s e m gene a es conse a i e (i.e., Liou ille) classical
dynamics and is gi en by
Lc=−p
m
∂
∂x+∂U(x)
∂x
∂
∂p.(2)
The second e m gene a es genuine quan um dynamics and is
gi en by
Lq=∞

n=1
(−1)n
(2n+1)!
¯h2n
4n
∂2n+1U(x)
∂x2n+1
∂2n+1
∂p2n+1.(3)
No e ha Lqis ze o o quad a ic po en ials (i.e., po en ials
wi h only linea and ha monic e ms). The hi d e m models
he p esence o noise and gene a es dissipa i e dynamics. Fo
le i a ed nanopa icles, i is con enien o conside [18]
Ln=γ1+p∂
∂p+¯h2
2x2

∂2
∂p2,(4)
whe e
=γkBT
¯h+1,(5)
kBis he Bol zmann cons an , and x=[¯h/(2m)]1/2is a
con enien leng h uni associa ed o he ze o-poin mo ion
luc ua ions o he quan um g ound s a e o a ha monic po-
en ial wi h equency . This sou ce o noise models a linea
coupling o a he mal ba h [19] o empe a u e T, wi h damp-
ing a e γ, and he p esence o a s ochas ic whi e- o ce e m
wi h displacemen noise a e gi en by 1.
The cen al poin o his pape is o use he Wigne unc ion
in a ime-dependen ame ha we call he Liou ille ame,
which is de ined as
˜
W(x,p, )≡e−Lc W(x,p, ).(6)
Since Lcis he gene a o o classical dynamics, one can use
he Liou ille heo em o w i e
˜
W(x,p, )=W(xc(x,p, ),pc(x,p, ), ),(7)
whe e xc(x,p, ) and pc(x,p, ) a e he solu ions o he clas-
sical equa ions o mo ion o poin pa icles mo ing in he
po en ial U(x) in he absence o noise wi h ini ial posi ion and
momen um gi en by xand p, espec i ely, namely, hey a e
solu ions o
∂xc(x,p, )
∂ =pc(x,p, )
m,
∂pc(x,p, )
∂ =−∂U(x)
∂xx=xc(x,p, )
,
(8)
wi h xc(x,p,0) =xand pc(x,p,0) =p. In he Liou ille
ame, he Wigne unc ion e ol es as
∂˜
W(x,p, )
∂ =e−Lc (Lq+Ln)eLc ˜
W(x,p, ).(9)
The Wigne unc ion in he Liou ille ame e ol es only
due o he p esence o quan um e ec s and/o noise, ha
is, ∂˜
W(x,p, )/∂ =0i Lq=Ln=0. The o iginal Wigne
unc ion W(x,p, ) can be ob ained om he Wigne unc ion
in he Liou ille ame ˜
W(x,p, )by
W(x,p, )=eLc ˜
W(x,p, )
=˜
W(xc(x,p,− ),pc(x,p,− ), ),
(10)
ha is, by using backwa d p opaga ion in ime o he classical
ajec o ies.
In he ollowing sec ion, we show ha nume ically sol ing
Eq. (9) on a ixed egula phase-space g id is highly e icien in
si ua ions in ol ing la ge expansions because i co esponds
o sol ing Eq. (1) on a ime-dependen , i egula phase-space
g id ha places g id poin s whe e hey a e mos c ucial. This
key idea is illus a ed in Fig. 1 o he example o a pa icle
e ol ing in a pu e qua ic po en ial, which we u he discuss
in Sec. IV.
III. NUMERICAL SIMULATION
IN THE LIOUVILLE FRAME
In his sec ion, we explain how o nume ically sol e he
ime e olu ion o he Wigne unc ion in he Liou ille ame,
namely, how o sol e Eq. (9). The i s s ep is o explic-
i ly calcula e he e ms in e−Lc (Lq+Ln)eLc . This allows
us o ob ain he explici o m o he pa ial de i a i e equa-
ion (PDE). As shown in he Appendix, one ob ains ha
Eq. (9) eads
∂˜
W(x,p, )
∂ =
n+m⩽NU

n,m=0
gnm(x,p, )∂n+m˜
W(x,p, )
∂xn∂pm.(11)
He e, NUis he smalles odd numbe such ha ∂nU(x)/∂xn=
0 o n⩾NU+2, which in u n de e mines ha Eq. (11)
is a PDE o o de NU. The ime-dependen scala unc ions
gnm(x,p, ) depend on he physical pa ame e s o he p oblem
(i.e., m,U(x), γ,T,1), bo h explici ly and implici ly h ough
he classical ajec o ies xc(x,p, ) and pc(x,p, ) and hei up
o NUo de de i a i es wi h espec o hei ini ial condi ion p.
013262-2
NUMERICAL SIMULATION OF LARGE-SCALE NONLINEAR … PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
FIG. 1. E olu ion in ime o he Wigne unc ion o a pa icle ini ially p epa ed in he g ound s a e o he ha monic po en ial Uh(x)and
e ol ing un il  =150 in he qua ic po en ial Uq(x)[seeEq.(16)] wi h η=102in he p esence o decohe ence wi h =10−5. Le panel
shows he ini ial s a e W(x,p,0) =˜
W(x,p,0), middle and igh panels show he e ol ed s a e in he o iginal ame W(x,p,150/)andin
he Liou ille ame ˜
W(x,p,150/), espec i ely. The black poin s, which appea as lines due o hei high densi y, depic a egula g id in he
Liou ille ame which we used o simula e he dynamics. The g id has 255 ×56 poin s wi h hx/x≈0.39 and hp/p≈0.16.
Thei de i a ion and explici exp essions o an up o qua ic
po en ial (NU=3) a e gi en in he Appendix.
The second s ep is o con e he PDE in Eq. (11)in oa
sys em o linea equa ions using he me hod o ini e di e -
ences. In he Liou ille ame, we use a uni o m ec angula
g id in xand pwi h sepa a ion be ween consecu i e g id
poin s gi en by hx>0 and hp>0 along each di ec ion, e-
spec i ely. The g id poin s a e gi en by (xi,pj)=(x0,p0)+
(ihx,jhp) o i=0,1,...,Nx−1 and j=0,1,...,Np−1.
He e, (x0,p0) is he bo om le poin o he g id which
con ains N=Nx×Nppoin s. The N alues o he Wigne
unc ion in he Liou ille ame ˜
W(x,p, ) e alua ed a he g id
poin s a e collec ed by he N-dimensional ec o ˜
W( ) whose
componen s, indexed by k=0,1,...,N−1, a e gi en by
˜
Wk=iNp+j( )=˜
W(xi,pj, ). Using a ini e di e ence me hod
(see he Appendix o u he de ails), one ob ains a sys em o
linea equa ions o his ec o gi en by
∂˜
W( )
∂ =D( )˜
W( ),(12)
whe e D( ) is a squa e N×Nma ix. Equa ion (12) can hen
be sol ed using
˜
W( + )=exp [D( ) ]˜
W( ),(13)
which is alid o a su icien ly small  (see he Appendix).
This nume ical me hod elies on de eloping a nume ically
e icien way o compu ing D( ), which equi es e alua ing
g(x,p, ) a he g id poin s. In u n, his equi es he al-
ues o xc(x,p, ) and pc(x,p, ) as well as he de i a i es
o xc(x,p,− ) and pc(x,p,− ) wi h espec o pa e e y
g id poin (xi,pj) and o all ins ances o ime conside ed
in he ini e di e ences app oach. Since an analy ical o -
mula o he classical ajec o ies and hei de i a i es wi h
espec o ini ial condi ions a e, in gene al, no a ailable o
nonha monic po en ials; hey need o be e icien ly e alua ed
nume ically. We do his by sol ing Eq. (8) and simila di e -
en ial equa ions ha can be de i ed o he de i a i es o he
classical ajec o ies wi h espec o ini ial condi ions using
a symplec ic me hod, which ensu es s abili y o e long in e-
g a ion imes [20]. Ano he impo an ool we use o imp o e
he e iciency o he me hod is o ela e he de i a i es o
he o wa d-p opaga ed ajec o ies wi h espec o he ini ial
condi ions wi h hose o he backwa d-p opaga ed ajec o ies
by making use o he p ope ies o he associa ed Jacobian
ma ices. We p o ide all he de ails o his nume ical me hod
in he Appendix, a me hod ha we ha e coded using C++,
Cy hon, and Py hon.
Fo mally, sol ing Eq. (11) in he ixed g id gi en by he
Nphase-space poin s (xi,pj) de ined abo e is equi alen o
sol ing Eq. (1) in a ime-dependen g id gi en by he Npoin s
(xi( ),pj( )) ≡(xc(xi,pj, ),pc(xi,pj, )), see Fig. 1.In his
ime-dependen g id, a ea u e which will be impo an o
ou la e discussion is he maximum phase-space g id den-
si y, namely, he minimal dis ance be ween wo phase-space
poin s. To quan i y his phase-space densi y, le us in o-
duce he ollowing dimensionless Jacobian ma ix o a gi en
phase-space poin , namely:
Jd(x,p, )≡⎛
⎝
∂xc(x,p, )
∂x
∂xc(x,p, )
∂p
p
x
∂pc(x,p, )
∂x
x
p
∂pc(x,p, )
∂p⎞
⎠,(14)
wi h p≡¯h/(2x). We hen de ine λ+
i,j( ) and λ−
i,j( )as
he la ges and smalles singula alues o he 2 ×2ma-
ix Jd(xi,pj, ). One can hen de ine a( )≡mini,jλ−
i,j( ),
namely, he minimum singula alue o e all he g id. In
his way, he dimensionless pa ame e 1/a( ) quan i ies he
maximum densi y o he phase-space ime-dependen g id. To
see his explici ly, conside a poin in phase space, w i en
wi hou dimensions as =(x/x,p/p) and a poin =
+(cos θ,sin θ) in i s close icini y (i.e., ||1). A e he
e olu ion go e ned by he classical ajec o y, he sepa a ion
be ween hese wo poin s can be exp essed, in linea o de in
013262-3
M. RODA-LLORDES e al. PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
,as
| c( , )− c( , )|
≈|Jd(x,p, )(cos θ,sin θ)T|.(15)
Acco ding o he singula alue decomposi ion o Jd(x,p, ),
he smalles singula alue o Jd(x,p, ) minimizes he dis-
ance Eq. (15) o e all possible di ec ions θ.
IV. EXAMPLE: QUARTIC POTENTIAL
Le us now apply he nume ical me hod p esen ed in his
pape o a pa icula example: quan um mechanics in a pu ely
qua ic po en ial [21]. A weak qua ic po en ial is capable
o p oducing la ge expansion and in e e ence inges on he
posi ion p obabili y dis ibu ion unc ion. This po en ial is
o in e es in he con ex o le i odynamics [1], especially
in expe imen s whe e nanopa icles e ol e in nonha monic
po en ials [4,6,22]. In addi ion, his example will show he
applicabili y o he nume ical me hod and illus a e ha sol -
ing ˜
W(x,p, ) in a cons an and egula g id is equi alen o
sol ing W(x,p, )inasma ime-dependen i egula phase-
space g id, see Fig. 1.
We conside a pa icle o mass m, whose s a e a =0,
namely, W(x,p,0) =˜
W(x,p,0), is gi en by he g ound s a e
o he ha monic po en ial Uh(x)=m2x2/2; see le panel
o Fig. 1. The posi ion and momen um s anda d de ia ion o
he ini ial s a e a e gi en by ˆx2=x=√¯h/(2m) and
ˆp2=p=¯h/(2x), espec i ely. A >0, he pa icle
e ol es in a pu ely qua ic po en ial U(x)=Uq(x), which we
pa ame ize as
Uq(x)=1
η4
¯h
4x
x4
.(16)
We conside he case o ic ionless noise (e.g., dynamics
in ul ahigh acuum [23]), namely, γ=0bu / > 0in
Eq. (4). The dimensionless pa ame e ηcha ac e izes he
s eng h o he qua ic po en ial. Since he ini ial kine ic en-
e gy o he s a e is ¯h/4, he u ning poin xdacco ding
o classical mechanics, de ined as ¯h/4=Uq(xd), is gi en
by xd/x=η. La ge phase-space expansions, namely, s a es
wi h spa ial delocaliza ion o de s o magni ude la ge han x
[24] will be hus gene a ed o weak po en ials, i.e., η1.
Le us i s analyze he e olu ion o he i s and second
phase-space momen s. Due o he alignmen o he qua ic po-
en ial wi h he ini ial s a e, he i s momen s emain cons an
and equal o ze o, namely, ˆx( )/x=ˆp( )/p=0. The
dynamics o he second momen s is shown in Fig. 2(a), whe e
we plo ˆx2( )/(ηx), ˆp2( )/(p), and {ˆx,ˆp}( )/(η¯h).
Using he η-scaled dimensionless imescale /η, hisplo
is o η1 con enien ly independen o η. The plo shows
ha he s a e expe iences ee dynamics du ing an ini ial
imescale gi en by 0 < /η 0.4, whe e ˆx2( )/(ηx)≈
{ˆx,ˆp}( )/(η¯h) g ows linea ly in ime and ˆp2( )/(p)
emains cons an and equal o one. A e his ini ial ime
in e al, he s a e s a s o expe ience he qua ic po en ial. In
pa icula , a /η ≈1.12, when {ˆx,ˆp}( )/(η¯h) is equal o
ze o, he s a e eaches a maximum alue o ˆx2( )/(ηx)
o he o de o 1, ha is, he s a e is spa ially delocalized o a
la ge leng h scale gi en by ηx[24]. This expansi e dynamics
FIG. 2. (a) Second-o de momen s o a pa icle ini ially p epa ed
in he g ound s a e o Uh(x) and e ol ing in Uq(x)wi hη=103
and and =2×10−8. The esul s o η=10,102, and 104in
hese no malized uni s is indis inguishable in he scale o he plo .
The e ical lines indica e he ins an whe e he ˆx2( )/(ηx)is
maximum and he ins an whe e {ˆx,ˆp}( )/(η¯h) eaches i s mos
nega i e alue, espec i ely. (b) G id densi y pa ame e a( )asa
unc ion o ime o he qua ic po en ial Uq(x) o wo di e en
alues o η.
ha gene a es a squeezed s a e is con enien ly accompanied
wi h an inc ease o he phase-space g id densi y. This can be
shown in Fig. 2(b), whe e we plo he η-scaled g id dis ance
ηa( ) as a unc ion o ime, showing ha he phase-space
densi y g ows as a unc ion o ime and is scaled wi h η.
Le us now s udy he e olu ion o he Wigne unc ion. In
Fig. 3,weshow ˜
W(x,p, ) (panel a) and W(x,p, ) (panel
b) o η=103and / =2×10−8a h ee ins ances o
ime: (i) =0, (ii)  /η ≈1.12 when ˆx2( )/(ηx)is
he la ges and he s a e gene a es an in e e ence pa e n in
he momen um p obabili y dis ibu ion, and (iii)  /η ≈1.56
when {ˆx,ˆp}( )/(2η¯h) eaches i s mos nega i e alue and
FIG. 3. Phase space ep esen a ion o he s a e o a pa icle ini-
ially p epa ed in he g ound s a e o Uh(x) and e ol ing inUq(x) wi h
η=103and =2×10−8. Resul s a di e en ins ances o ime,
namely, he ini ial s a e, he ime whe e he a iance in posi ion is
maximum, and he ime when he co a iance is minimum, see Fig. 2.
(a) ˜
W(x,p, ). In his case, we used a g id wi h 2048 ×256 poin s
wi h hx/x≈0.24 and hp/p≈0.06. (b) W(x,p, ).
013262-4
NUMERICAL SIMULATION OF LARGE-SCALE NONLINEAR … PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
FIG. 4. (a) P obabili y dis ibu ion in posi ion P(x) a ime
/η =1.56 o a pa icle ini ially p epa ed in he g ound s a e o
Uh(x) and e ol ing in Uq(x) wi h η=103and / =2×10−8.This
ime co esponds o he momen whe e he co a iance is minimum,
see Fig. 2. The inse shows how he sepa a ion be ween he i s
wo peaks x scales as a unc ion o η. (b) Visibili y o he second
la ges maximum o P(x) a he ime speci ied abo e, as a unc ion
o displacemen noise a e and o di e en alues o η.
he s a e exhibi s an in e e ence pa e n in he posi ion p ob-
abili y dis ibu ion, see Fig. 4(a).InFig.2(a), he ins ances o
ime (ii) and (iii) a e indica ed wi h a e ical dashed line. We
emphasize ha W(x,p, ) [Fig. 2(b)] is ob ained by simply
using Eq. (10) a e ha ing nume ically ob ained ˜
W(x,p, )
[Fig. 2(a)] wi h he me hod p esen ed in his pape . Compa ing
he xaxes o Figs. 3(a) and 3(b), one can see how W(x,p, )
expands signi ican ly mo e han ˜
W(x,p, ). As shown in
Fig. 1, he egula g id poin s used o ep esen ˜
W(x,p, )in
he Liou ille ame a e e icien ly dis ibu ed in he o iginal
ame o p ope ly desc ibe W(x,p, ). The esul s in Fig. 3(a)
a e ob ained in a ixed g id o a phase-space leng h scale
gi en by (hx/x)2+(hp/p)2≈0.25. The leng h scale in
he ime-dependen g id, namely, in Fig. 3(b) is educed by a
ac o o a( ), which as one can see in Fig. 2(b), eaches alues
below 10−3, way below he phase-space Planck scale [15].
This means ha o ma ch he accu acy le el o ou me hod,
using a egula g id in he o iginal ame would need abou
103 imes mo e g id poin s in each di ec ion. This illus a es
he ad an age p o ided by ou app oach.
Finally, le us discuss he impac o noise by illus a -
ing how i a ec s he isibili y o he in e e ence pa e n
in posi ion a he ime /η ≈1.56. In Fig. 4(a),weplo
he p obabili y dis ibu ion P(x)≡∞
−∞ dpW(x,p, )a his
pa icula ins ance o ime o η=103and / =2×10−8.
We de ine x as he dis ance be ween he la ges in e e ence
peak and i s neighbo ing peak. Fo he pa ame e s in Fig. 4,we
ob ain x /x≈21.2. As shown in he inse o Fig. 4(a), he
scaling o his dis ance wi h ηis gi en by x /x≈2.11η1/3.
The isibili y o his in e e ence pa e n, de ined as (Pmax −
Pmin )/(Pmax +Pmin ), whe e Pmax and Pmin a e he alues o
P(x) a he la ges maximum and i s neighbo ing minimum,
espec i ely, is a dec easing a unc ion o / as we show in
Fig. 4(b). As expec ed [6,24], he impac o in he isibili y
scales oughly as η2.
The s udy o quan um dynamics in a nonha monic po-
en ial in he p esence o noise, which we ha e pe o med
using he nume ical me hod p esen ed in his pape , is ele an
o cu en e o s o p epa e la gely delocalized mac oscopic
quan um s a es o la ge masses [4,6,25–32]. We ema k ha ,
easibili ywise, pu ely qua ic po en ials a e no ideal since he
imescale needed o gene a e he in e e ence pa e n shown
in Eq. (4), ha is, /η ≈1.56, is o η1 much la ge
han he a e age collision ime wi h a single gas molecule
a ul ahigh acuum [18]. This is one o he main easons
mo i a ing ou ecen p oposal [6], which is also analyzed
wi h he nume ical me hod p esen ed in his pape , whe e a
double-well po en ial such ha he in e ed ha monic e m
exponen ially speeds up he dynamics [5,32]isused.
V. CONCLUSIONS
In his pape , we ha e p esen ed a me hod o simula e
nonlinea open quan um dynamics o po en ials in which
he quan um s a e expands se e al o de s o magni ude in
phase space while exhibi ing ele an ea u es a e y small
sub-Planck scales [15]. This egime is o pa icula in e es
o designing, op imizing, and unde s anding p o ocols ha
gene a e mac oscopic quan um s a es by le ing a massi e
pa icle e ol e in a nonha monic po en ial [3,6,33]. We ha e
demons a ed he powe o his me hod using he dynamics
o an ini ially highly localized s a e in a qua ic po en ial. We
ha e shown how in his po en ial he s a e posi ion a iance
g ows by se e al o de s o magni ude, and ye i s Wigne
unc ion exhibi s nega i e ea u es on a scale below he ini ial
ze o-poin luc ua ions. P ope ly desc ibing such small scales
using a egula g id in he o iginal ame would equi e an
imp ac icable amoun o poin s, a challenge ha we o e come
by he in oduc ion o he Liou ille ame.
This me hod should be applicable o a b oade class o
quan um mechanical p oblems. In p inciple, any po en ial
U(x) can be conside ed, howe e , he numbe o de i a i es
conside ed in Eq. (9) mus be ini e o allow o a nume ical
e alua ion. In oducing a cu o o he o de o he po en-
ial in U(x) should yield accu a e esul s. In addi ion, o he
ypes o noise and decohe ence beyond he ones conside ed
in his pape (e.g., s ochas ic o ce-g adien ) can also be
inco po a ed. Fu he mo e, while we ha e conside ed bo h
ime-independen po en ials and decohe ence a es, he nu-
me ical me hod is inhe en ly ime dependen , see Eq. (13),
which means ha ime dependence could be in oduced wi h
he co esponding modi ica ions. An ad an ageous ea u e o
s udying quan um mechanics wi h he Wigne unc ion is
ha he classical limi can be easily aken, namely, aking
¯h=0 such ha Lq=0inEq.(1). In his classical limi ,
he Wigne unc ion in he Liou ille ame is only d i en by
dissipa i e dynamics. We ha e ocused on a one-dimensional
p oblem, bu he me hod can be gene alized o highe spa ial
dimensions. Finally, while we ha e ocused on he mo ion
013262-5

M. RODA-LLORDES e al. PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
o a pa icle in wide nonha monic po en ials, ou me hod is
gene al and could also be applied o s udy simila nonlinea
dynamics o a bosonic mode in o he pla o ms. A ield whe e
i could be ele an is quan um in o ma ion p ocessing wi h
ci cui elec odynamics, see Re . [34] and e e ences he ein.
In conclusion, he me hod p esen ed in his pape elies
on a c ucial elemen : he desc ip ion o Wigne unc ion dy-
namics in he Liou ille ame Eq. (9). We emphasize ha his
ame p o es o be highly aluable no only in p ac ical e ms
bu also om a concep ual s andpoin , as i clea ly un eils
he impac o quan um physics in he mechanical mo ion o
a pa icle.
ACKNOWLEDGMENTS
We would like o hank C. Dellago, L. Einkemme , D.
Giannand ea, M. Inne bichle , T. Weiss, and he Q-X eme
syne gy g oup o help ul discussions. This esea ch has been
suppo ed by he Eu opean Resea ch Council (ERC) unde
G an Ag eemen No. 951234 (Q-X eme ERC-2020-SyG)
and by he Eu opean Union’s Ho izon 2020 esea ch and inno-
a ion p og am unde G an Ag eemen No. 863132 (IQLe ).
P.T.G. was pa ially suppo ed by he Founda ion o Polish
Science (FNP).
APPENDIX: DETAILS ON THE NUMERICAL METHOD
In his Appendix, we de ail all he s eps used o sol e
Eq. (9) nume ically.
1. Explici exp ession o he PDE
The i s s ep is o ob ain an explici exp ession o
e−Lc (Lq+Ln)eLc . To ind how any ope a o O ans o ms
unde e−Lc Oe
Lc , one can apply he ans o med ope a o
o an a bi a y unc ion (x,p) and iden i y which ope a o
p oduces he same esul . I is use ul o ecall ha by i ue o
he Liou ille heo em, we know how e±Lc ac s on an a bi a y
unc ion (x,p), namely,
e±Lc (x,p)= (xc(x,p,∓ ),pc(x,p,∓ )).(A1)
In ou case, he ope a o s ha appea in Lq+Lna e x,p,
and de i a i es wi h espec o p.Fo xand p, making use o
Eq. (A1), one inds ha hese ope a o s ans o m acco ding
o
e−Lc xe
Lc =xc(x,p, ) and e−Lc pe
Lc =pc(x,p, ).
(A2)
Simila ly, one inds ha he de i a i e wi h espec o p ans-
o ms acco ding o he chain ule as
e−Lc ∂
∂peLc =∂p¯xc
∂
∂x+∂p¯pc
∂
∂p,(A3)
whe e we in oduce he ollowing sho hand no a ion:
∂n
p¯xc=∂nxc(x,p,− )
∂pnx=xc(x,p, )
p=pc(x,p, )
and ∂n
p¯pc=∂npc(x,p,− )
∂pnx=xc(x,p, )
p=pc(x,p, )
.(A4)
No e ha hese a e scala unc ions o x,p, and , and hey a e he de i a i es wi h espec o ini ial condi ions o he classical
ajec o ies s a ing om he poin (xc(x,p, ),pc(x,p, )) p opaga ed backwa ds in ime o a ime . Explici ly, o n=1 hey
co espond o he ollowing limi :
∂p¯xc=lim
ε→0
xc(xc(x,p, ),pc(x,p, )+ε, − )−xc(xc(x,p, ),pc(x,p, ),− )
ε=lim
ε→0
xc(xc(x,p, ),pc(x,p, )+ε, − )−x
ε.
(A5)
Fo highe o de de i a i es, one inds exp essions co esponding o mul iple applica ions o he chain ule, namely, o he
second-o de de i a i e wi h espec o p, one has
e−Lc ∂2
∂p2eLc =∂2
p¯xc
∂
∂x+∂2
p¯pc
∂
∂p+(∂p¯xc)2∂2
∂x2+(∂p¯pc)2∂2
∂p2+2∂p¯xc∂p¯pc
∂2
∂x∂p,(A6)
whe eas, o he hi d-o de de i a i e, one has
e−Lc ∂3
∂p3eLc =∂3
p¯xc
∂
∂x+∂3
p¯pc
∂
∂p+3∂p¯xc∂2
p¯xc
∂2
∂x2+∂p¯pc∂2
p¯pc
∂2
∂p2+∂2
p¯xc∂p¯pc+∂p¯xc∂2
p¯pc∂2
∂x∂p
+(∂p¯xc)3∂3
∂x3+(∂p¯pc)3∂3
∂p3+3(∂p¯xc)2∂p¯pc
∂3
∂x2∂p+∂p¯xc(∂p¯pc)2∂3
∂x∂p2.(A7)
We only conside po en ials U(x) o which i h and highe
o de de i a i es anish, and he e o e only de i a i es up o
hi d o de will appea in Eq. (9). Fo po en ials whe e highe
o de de i a i es a e ele an , one could ex end ou app oach
o include hem. Subs i u ing Eqs. (A2), (A3), (A6), and (A7)
in o Eq. (9) yields he explici equa ion ha we need o sol e
nume ically. I has he ollowing o m:
∂˜
W(x,p, )
∂ =
n+m⩽3

n,m=0
gnm(x,p, )∂n+m˜
W(x,p, )
∂xn∂pm,(A8)
013262-6
NUMERICAL SIMULATION OF LARGE-SCALE NONLINEAR … PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
whe e he explici exp essions o he coe icien s a e gi en by
g00(x,p, )=γ, (A9)
g10(x,p, )=γpc∂p¯xc+¯h2
2x2

∂2
p¯xc−¯h2
12U(3)(xc)∂3
p¯xc,
(A10)
g01(x,p, )=γpc(x,p, )∂p¯pc+¯h2
2x2

∂2
p¯pc
−¯h2
12U(3)(xc)∂3
p¯pc,(A11)
g20(x,p, )=¯h2
2x2

(∂p¯xc)2−¯h2
4U(3)(xc)∂p¯xc∂2
p¯xc,(A12)
g02(x,p, )=¯h2
2x2

(∂p¯pc)2−¯h2
4U(3)(xc)∂p¯pc∂2
p¯pc,(A13)
g11(x,p, )=¯h2
2x2

∂p¯xc∂p¯pc−¯h2
4U(3)(xc)
×∂p¯xc∂2
p¯pc+∂p¯pc∂2
p¯xc,(A14)
g30(x,p, )=−¯h2
12U(3)(xc)(∂p¯xc)3,(A15)
g03(x,p, )=−¯h2
12U(3)(xc)(∂p¯pc)3,(A16)
g21(x,p, )=−¯h2
4U(3)(xc)(∂p¯xc)2∂p¯pc,(A17)
g12(x,p, )=−¯h2
4U(3)(xc)(∂p¯pc)2∂p¯xc.(A18)
To simpli y no a ion, he e and he ea e we use U(i)(x) o he
i h de i a i e o Ue alua ed a x. Also, no e ha we use xc
and pcas a sho hand o xc(x,p, ) and pc(x,p, ) o simpli y
he exp essions, bu hey s ill depend on x,p, and .
2. Disc e iza ion o he PDE
Now ha we ha e an explici exp ession o he equa-
ion we need o sol e, we need o disc e ize i o allow o
nume ical simula ion. To do so, we desc ibe ˜
Win a egula
g id which con ains N=Nx×Nppoin s which we deno e by
(xi,pi). Then, we deno e he alues o ˜
Win each o hese g id
poin s by ˜
Wi,j=˜
W(xi,pj). Nex , we exp ess he de i a i es
wi h espec o xand pin Eq. (A8) in e ms o ini e di e ence
schemes. In pa icula , we use a second-o de cen e ed ini e
di e ence scheme, which we lis below o he i s -, second-,
and hi d-o de de i a i es. Fi s , o he i s -o de de i a i es
hey ead
∂˜
Wi,j
∂x=
˜
Wi+1,j−˜
Wi−1,j
2hx
,(A19)
∂˜
Wi,j
∂p=
˜
Wi,j+1−˜
Wi,j−1
2hp
.(A20)
Nex , o he second-o de de i a i es, one has
∂2˜
Wi,j
∂x2=
˜
Wi+1,j+˜
Wi−1,j−2˜
Wi,j
h2
x
,(A21)
∂2˜
Wi,j
∂p2=
˜
Wi,j+1+˜
Wi,j−1−2˜
Wi,j
h2
p
,(A22)
∂2˜
Wi,j
∂x∂p=
˜
Wi+1,j+1+˜
Wi−1,j−1−˜
Wi−1,j+1−˜
Wi+1,j−1
4hxhp
.
(A23)
Finally, he exp essions o he hi d-o de de i a i es a e
gi en by
∂3˜
Wi,j
∂x3=
˜
Wi+2,j−2˜
Wi+1,j+2˜
Wi−1,j−˜
Wi−2,j
2h3
x
,(A24)
∂3˜
Wi,j
∂p3=
˜
Wi,j+2−2˜
Wi,j+1+2˜
Wi,j−1−˜
Wi,j−2
2h3
p
,(A25)
∂3˜
Wi,j
∂x2∂p=
˜
Wi+1,j+1+˜
Wi−1,j+1−2˜
Wi,j+1+2˜
Wi,j−1−˜
Wi+1,j−1−˜
Wi−1,j−1
2h2
xhp
,(A26)
∂3˜
Wi,j
∂x∂p2=
˜
Wi+1,j+1+˜
Wi+1,j−1−2˜
Wi+1,j+2˜
Wi−1,j−˜
Wi−1,j+1−˜
Wi−1,j−1
2hxh2
p
.(A27)
A e subs i u ing all he de i a i es in Eq. (A8) by hei ini e
di e ence e sions [see Eqs. (A19)–(A27)], he igh -hand
side o he equa ion is gi en by a linea combina ion o ˜
Wi,j
wi h di e en indices i,j. Explici ly, one has
∂˜
Wi,j( )
∂ =
α
D(i,j),(α)˜
Wα( ),(A28)
whe e α uns o e he ollowing 13 indices: (i,j),(i±
1,j),(i,j±1),(i±1,j±1),(i±1,j∓1),(i±2,j) and
(i,j±2). Collec ing he alues o ˜
Wi,jin he N-dimensional
ec o ˜
W( ) wi h componen s ˜
Wk=iNp+j( )=˜
Wi,j( ) indexed
by k=0,1,...,N−1 allows us o w i e Eq. (A28)as
∂˜
W( )
∂ =D( )˜
W( ),(A29)
whe e D( )isaN×Nma ix. F om his equa ion, one
can de i e an exp ession o p opaga e he solu ion in ime
013262-7
M. RODA-LLORDES e al. PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
gi en by
˜
W( + )=exp  +
D( )d ˜
W( )
≈exp [D( ) ]˜
W( ),(A30)
whe e he app oxima ion assumes ha  is small enough
such ha D( ) a ies slowly enough be ween and + .
The en ies o he D( ) ma ix can be ound by inspec-
ion a e eplacing he de i a i es in Eq. (A8) by hei ini e
di e ence e sions [see Eqs. (A19)–(A27)]. Fo ins ance, he
ma ix en y co esponding o he index (i,j),(i+1,j) eads
D(i,j),(i+1,j)( )=g10(xi,pj, )
2hx+g20(xi,pj, )
h2
x
−g30(xi,pj, )
h3
x−g21(xi,pj, )
hxh2
p
.(A31)
No ice ha each ow o D( ) will only ha e 13 en ies di e -
en om ze o, which means ha D( ) will be spa se. This is
due o he ac ha ini e di e ences only ela e poin s wi h
up o second-o de neighbo s. One could ha e chosen highe -
o de ini e di e ences, in which case he e would me mo e
nonze o en ies in each ow o D( ). Howe e , we ound ha
inc easing he ini e di e ences om second o ou h o de
didn’ yield any signi ican imp o emen in he accu acy o
ou solu ion. Finally, no e ha o ully de ine D( ), one needs
o speci y he bounda y condi ions. We use pe iodic bounda y
condi ions since hey p o ide a mo e s able simula ion han
ze o- alue bounda y condi ions. In pa icula , we iden i y he
igh and op edges o he g id wi h he le and bo om edges,
espec i ely. Explici ly, we iden i y i=Nxwi h i=0, and
j=Npwi h j=0.
3. E icien compu a ion o he Dma ix
As one can see in Eq. (A31), ob aining he nume ical alue
o he di e en en ies o he D( ) ma ix equi es e alu-
a ing all gmn(xi,pj, )[seeEqs.(A9)–(A18)] in each poin
o he g id. In u n, his equi es he alues o xc(xi,pj, )
and pc(xi,pj, ) as well as he de i a i es ∂n
p¯xcand ∂n
p¯pc[see
Eq. (A4)] up o n=3 a e e y g id poin (xi,pj) and o all
ins ances o ime conside ed in he ini e di e ences app oach.
Since an analy ical o mula o he classical ajec o ies is gen-
e ally no a ailable o nonha monic po en ials we e alua e
hem nume ically.
We ob ain xc(xi,pj, ) and pc(xi,pj, ) by p opaga ing in
ime he classical equa ions o mo ion Eq. (8) wi h each g id
poin (xi,pj) as ini ial condi ion. To ensu e s abili y o e long
in eg a ion imes we use a symplec ic me hod [20]. In pa icu-
la , we use he ou h-o de me hod desc ibed in Re . [35]. To
ob ain he de i a i es o he in e se mapping ∂n
p¯xcand ∂n
p¯pc,
we use an app oach consis ing o wo s eps. Fi s , we compu e
he de i a i es o he di ec mapping as solu ions o di e en-
ial equa ions, which allows us o bene i om he p ope ies
o he symplec ic me hod used abo e. Second, we use hese
alues o compu e ∂n
p¯xcand ∂n
p¯pc h ough he ela ion be ween
he di ec and in e se mapping. Using hese s eps is mo e
e icien han a di ec nume ical e alua ion o hese de i a i es
in e ms o limi s such as he one shown in Eq. (A5). In he
ollowing, we desc ibe hese wo s eps in de ail.
By aking de i a i es wi h espec o xand pin Eq. (8),
one can ob ain he equa ion o mo ion o he de i a i es we
need. No e ha we use xcand pcas sho hand o xc(x,p, )
and pc(x,p, ), espec i ely. Speci ically, aking he de i a i e
wi h espec o xon Eq. (8) yields he di e en ial equa ions o
∂xxc(x,p, ) and ∂xpc(x,p, ):
∂
∂xc
∂x=1
m
∂pc
∂x,∂
∂pc
∂x=−U(2)(xc)∂xc
∂x.(A32)
The ini ial condi ions a e gi en by ∂xxc(x,p,0) =1 and
∂xpc(x,p,0) =0. They s em om he ac ha , a ime =
0, xc(x,p,0) =xand xc(x,p,0) =p. Simila ly, aking he
de i a i e wi h espec o pyields a simila equa ion o
∂pxc(x,p, ) and ∂ppc(x,p, ),
∂
∂xc
∂p=1
m
∂pc
∂p,
∂
∂pc
∂p=−U(2)(xc)∂xc
∂p,
(A33)
wi h ini ial condi ions ∂pxc(x,p,0) =0 and
∂ppc(x,p,0) =1. By aking mo e de i a i es, one can ob ain
equa ions o he highe o de de i a i es. The second-o de
de i a i es wi h espec o ini ial condi ions ul ill
∂
∂2xc
∂x2=1
m
∂2pc
∂x2,
∂
∂2pc
∂x2=−U(3)(xc)∂xc
∂x2
−U(2)(xc)∂2xc
∂x2,
(A34)
∂
∂2xc
∂p2=1
m
∂2pc
∂p2,
∂
∂2pc
∂p2=−U(3)(xc)∂xc
∂p2
−U(2)(xc)∂2xc
∂p2,
(A35)
∂
∂2xc
∂x∂p=1
m
∂2pc
∂x∂p,
∂
∂2pc
∂x∂p=−U(3)(xc)∂xc
∂x
∂xc
∂p−U(2)(xc)∂2xc
∂x∂p,
(A36)
wi h all he ini ial condi ions being ze o. The hi d-o de
de i a i es ul ill
∂
∂3xc
∂x3=1
m
∂3pc
∂x3,
∂
∂3pc
∂x3=−U(4)(xc)∂xc
∂x3
−3U(3)(xc)∂xc
∂x
∂2xc
∂x2
−U(2)(xc)∂3xc
∂x3,(A37)
∂
∂3xc
∂p3=1
m
∂3pc
∂p3,
∂
∂3pc
∂p3=−U(4)(xc)∂xc
∂p3
−3U(3)(xc)∂xc
∂p
∂2xc
∂p2
−U(2)(xc)∂3xc
∂p3,(A38)
013262-8
NUMERICAL SIMULATION OF LARGE-SCALE NONLINEAR … PHYSICAL REVIEW RESEARCH 6, 013262 (2024)
∂
∂3xc
∂x2∂p=1
m
∂3pc
∂x2∂p,
∂
∂3pc
∂x2∂p=−U(4)(xc)∂xc
∂x2∂xc
∂p−U(3)(xc)∂2xc
∂x2
∂xc
∂p
−2U(3)(xc)∂xc
∂x
∂2xc
∂x∂p−U(2)(xc)∂3xc
∂x2∂p,(A39)
∂
∂3xc
∂x∂p2=1
m
∂3pc
∂x∂p2,
∂
∂3pc
∂x∂p2=−U(4)(xc)∂xc
∂x∂xc
∂p2
−U(3)(xc)∂xc
∂x
∂2xc
∂p2
−2U(3)(xc)∂xc
∂p
∂2xc
∂x∂p−U(2)(xc)∂3xc
∂x∂p2,(A40)
again wi h all he ini ial condi ions being ze o.
No e ha xcand pcappea explici ly in all equa ions. Sim-
ila ly, ∂xxc,∂pxc,∂xpc, and ∂ppcappea in he equa ions o
he second- and hi d-o de de i a i es, and he second-o de
de i a i es appea in he equa ions o he hi d-o de de i a-
i es. This means ha o sol e he equa ions o highe o de
de i a i es, he alues o all he lowe de i a i es a e needed
as an inpu . E en mo e, no only he alues a each ime
being conside ed a e needed, bu also he alues a he ou
in e media e ime s eps in he ou h-o de me hod [35] ha
we use. To be memo y e icien , we do no use a sepa a e
sol e o each equa ion bu a he a single sol e o all equa-
ions ha co ec ly use all he p e iously compu ed alues in
he igh sequence.
Finally, we need o ela e hese de i a i es o he de i a-
i es o he in e se map ∂n
p¯xcand ∂n
p¯pc. Fo he i s -o de
de i a i es, he key obse a ion is ha he Jacobian ma ix o
he map
J(x,p, )≡⎛
⎝
∂xc(x,p, )
∂x
∂xc(x,p, )
∂p
∂pc(x,p, )
∂x
∂pc(x,p, )
∂p,⎞
⎠(A41)
is by cons uc ion he in e se o he Jacobian ma ix o he
in e se map:
˜
J(x,p, )=∂x¯xc∂p¯xc
∂p¯xc∂p¯pc.(A42)
Using his ac , we compu e J(xi,pj, ) o each poin in he
g id a each ime s ep, and hen ob ain ˜
Jby in e ing he ma-
ix. Explici ly, we use he ollowing o mula: ˜
J(xi,pj, )=
J−1(xi,pj, ). One can show ha he de e minan o bo h ˜
J
and Jis cons an and equal o one and, he e o e, compu ing
his in e se is s aigh o wa d. Simila ela ionships exis o
highe o de de i a i es, which we de i e below.
To simpli y he exp essions in he ollowing, we de-
ine he ec o =(x,p) and he ec o unc ion c( , )=
(xc( , ),pc( , )). Finally, we de ine a new se o a iables
˜ =(˜x,˜p) which a e ela ed o h ough he classical ajec-
o ies as
= c(˜ , ) o , equi alen ly, ˜ = c( ,− ).(A43)
Using his no a ion, we can exp ess he Jacobian ma ices
discussed abo e as
Ji
j=∂ i
∂˜ j
and ˜
Ji
j=∂˜ i
∂ j
,(A44)
and hei ela ionship o being he in e se o each o he as
kJi
k˜
Jk
j=δij. Now, o de i e a ela ion o he second-o de
de i a i es, we s a by de ining he Hessian enso and in e se
Hessian enso , espec i ely, as
Hi
jk =∂2 i
∂˜ j∂˜ k
and ˜
Hi
jk =∂2˜ i
∂ j∂ k
,(A45)
whe e i,j,kcan be ei he 1 o 2. Nex , we expand he ollow-
ing exp ession using he chain ule:
0=∂2 i
∂ j∂ k=∂
∂ k
l
∂ i
∂˜ l
∂˜ l
∂ j
=
l
∂ i
∂˜ l
∂2˜ l
∂ j∂ k+
l,m
∂2 i
∂˜ l∂˜ m
∂˜ l
∂ j
∂˜ m
∂ k
.(A46)
Then, using he p ope ies o he Jacobian ma ices, we can
ew i e he exp ession abo e as
0=
l
Ji
l˜
Hl
j,k+
l,m
Hi
l,m˜
Jl
j˜
Jm
k
o ˜
Hi
j,k=−
n,l,m
Hn
l,m˜
Ji
n˜
Jl
j˜
Jm
k.(A47)
One can hen use his exp ession o ob ain he alues o ˜
H
in e ms o H(which we compu e by sol ing he di e en ial
equa ions desc ibed abo e) and he alues o ˜
J ha we al eady
compu ed.
Fo he hi d-o de de i a i es one can p oceed in a simila
ashion. One de ines he enso s
˜
Ti
jkα=∂3˜ i
∂ j∂ k∂ α
and Ti
jkα=∂3 i
∂˜ j∂˜ k∂˜ α
(A48)
and akes ye ano he de i a i e wi h espec o αin
Eq. (A46). Then, p oceeding in a simila way, one inally
a i es a he exp ession
˜
Ti
j,k,α =− 
n,l,m,β
Tn
l,m,β ˜
Ji
n˜
Jl
j˜
Jm
k˜
Jβ
α
−
n,l,m
Hn
l,m˜
Ji
n˜
Hl
j,k˜
Jm
α+˜
Hl
j,α ˜
Jm
k+˜
Hl
k,α ˜
Jm
j.(A49)
In summa y, ou nume ical app oach o sol e Eq. (11)
consis s o he ollowing s eps o each ime s ep  . Fi s ,
p opaga e in ime he classical ajec o ies, and i s de i a i es
wi h espec o ini ial condi ions, o each poin in he g id.
Second, use hese de i a i es o compu e he co esponding
de i a i es o he in e se map. Thi d, use all hese newly com-
pu ed alues o gene a e he ma ix D( ). Finally, use Eq. (13)
o compu e ˜
Wa he new ime s ep in e ms o he alues a
he p e ious ime s ep. Repea ing his p ocedu e allows us
o p opaga e ˜
Win ime. We implemen ed all hese s eps by
de eloping ou own simula ion code in C++, Cy hon, and
Py hon.
013262-9