Numerical simulation of large-scale nonlinear open quantum mechanics

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