Many theories of quantum gravity can be understood as imposing a minimum length scale the signatures of which can potentially be seen in precise table top experiments. In this work, we inspect the capacity for correlated many-body systems to probe non-classicalities of spacetime through modifications of the commutation relations. We find an analytic derivation of the dynamics for a single mode light field interacting with a single mechanical oscillator and with coupled oscillators to first order corrections to the commutation relations. Our solution is valid for any coupling function as we work out the full Magnus expansion. We numerically show that it is possible to have superquadratic scaling of a nonstandard phase term, arising from the modification to the commutation relations, with coupled mechanical oscillators.

## I. INTRODUCTION

Establishing a quantum theory of gravity is one of the greatest challenges in contemporary physics. With early indications originating from gravitational instabilities arising in ultra-high-energy scattering in string theory,^{1–3} it is now widely thought in many different approaches to quantum gravity that there should be a minimum length scale,^{4,5} typically taken to be on the order of the Planck length $lp=\u210fG/c3$, where $\u210f$ is the reduced Planck constant, *G* is the gravitational constant, and *c* is the speed of light in vacuum. This minimum length scale can be arrived at from general heuristic arguments,^{6} such as the fact that a probe signal that could resolve such lengths would have an associated stress-energy distribution that collapses it into a black hole. Because the standard Heisenberg uncertainty relations do not prevent arbitrary fine squeezing of a state in position space, at the expense of high uncertainty in momentum space, the existence of a minimum length scale comes with a generalized uncertainty principle (GUP) $\Delta x\Delta p\u2265(\u210f/2)(1+\beta \Delta p2)$, where $\beta =\beta 0/(Mpc)2$ is what we shall call the GUP correction,^{7} and *M _{p}* is the Planck mass.

^{8,9}In turn this can be viewed, using Robertson's inequalities,

^{10}as resulting from a typical deformation of canonical commutation relations to the nonstandard form $[xi,pj]=i\u210f(1+\beta pj2)\delta ij$.

^{11}It has been pointed out recently

^{12}that this algebraic structure of the GUP may not be the most general, but the problem we address in this work also concerns more general forms. Regardless of the choice of the nonstandard commutator, there has yet to be empirical evidence suggesting any corrections from standard commutation relations, as current observations

^{13,14}do not rule out values of

*β*

_{0}much below 10

^{6}, while it is assumed

*β*

_{0}should be on the order of unit for Planck scale corrections.

^{15}Although non-classical features of spacetime and gravitation may be censored by some Penrose type objective collapse process due to gravity itself,

^{16–18}in recent years much attention has been given to proposals seeking to uncover some non-classicalities of spacetime while circumventing the requirement to work at the Planck energy level.

^{19–23}By relinquishing high-energy probes for directly observable Planckian effects, those proposals rely instead on finely controlled low-energy experiments that would allow us to indirectly detect signatures of deviations from the standard theory.

Here, we study a toy model inspired from a proposal that we term “optomechanical protocol”^{22} as well as a reaction to all protocols involving massive systems to test deformed commutators.^{24} The optomechanical protocol suggests that one can gain information on the type of evolution undergone by a mechanical oscillator by inspecting the non-linear phase that the light picks up after having interacted with the oscillator. A sequence of light pulses drives the mechanical oscillator around a close phase space trajectory, returning the oscillator to its initial state (in the absence of commutator deformations), while imparting a (mechanical state-independent) non-linear optical phase. This phase is of interest as it may carry nonstandard signatures and, as we shall show from first principles, is quartic in the photon number if GUP corrections exist, and is quadratic otherwise.

A major theoretical challenge to these types of proposals, regardless of experimental difficulties, is that one may, in principle, not be able to use the center of mass coordinates of any macroscopic systems, which are collections of elementary particles, to probe deformations of the commutation relations. Following Refs. 22 and 24, if one considers the center of mass coordinates of a system of *N* equally massive subsystems $X=(1/N)\u2211k=1Nxk$ and $P=\u2211k=1Npk$ with a deformed commutator $[xi,pj]=i\u210f(1+\beta pj2)\delta ij$, then

where we have expanded and simplified the commutator $(1/N)\u2211i,j[xi,pj]$. In the quasi-rigid limit ($\u2200i,\u2009pi\u2248P/N$), the variance of the momentum distribution over the subsystems vanishes, so that $[X,P]\u2248i\u210f(1+(\beta /N2)P2)$, and thus, deformations to the center of mass commutator are suppressed by the (potentially macroscopic) number of constituent particles. This leads to an important question: Are many-body macroscopic systems relevant probes for the Planck scale behavior?

Although the scaling of the suppression is unknown in general,^{25} from Eq. (1) one can expect a scaling between $1/N$ (for uncorrelated systems) and $1/N2$ (for rigid systems). However, in specially correlated fine-tuned systems, it may be that some Planck scale corrections could become amplified: this would be the case if some signature of the GUP correction is found to have a superquadratic scaling with the number of elementary constituents. In our work, this nonstandard signature arising from the interaction between a single mode optical field with a many-body system will be taken to be a phase term that is quartic in the photon number and is independent of the state of the many-body mechanical system. In order to observe superquadratic scaling, one must first be able to fully characterize the dynamics of many-body systems in a sufficiently general manner to allow for a broad type of interaction patterns and couplings between the subsystems. This accurate description is what we undertake in this work, and what we use to show numerical evidence of superquadratic scaling.

This paper is structured as follows. In Sec. II, we present the basic situation on which we build the rest of our work; we go through classical calculations in the case of a single particle interacting with a single mode optical field and show how the dynamics is solved with an order 2 Magnus expansion. In Sec. III, we modify the commutation relations and completely solve the full time-dependent dynamics for arbitrary interaction functions, to first order in the GUP corrections. We observe that the order 5 Magnus expansion gives the exact propagator and includes a phase term that is independent of the state of the mechanical oscillator, with quartic dependence on the photon number. In Sec. IV, we impose a specific form for the interaction function and show how to recover the phase term from first principles. In Sec. V, we extend our analysis to many-body systems and lay out a general derivation of the fifth Magnus term. In Sec. VI, we apply our general formula to a many-body system interacting with pulsed light and demonstrate a numerical example of superquadratic scaling behavior when the system is a chain of coupled harmonic oscillators.

## II. SINGLE PARTICLE TOY MODEL

Let us consider a single particle of mass *m* that is harmonically trapped with angular frequency Ω, interacting with a single mode optical field with an annihilation operator *a* and the angular frequency Ω_{F}. In standard quantum mechanics, the optomechanical Hamiltonian reads

where *x* is the position of the particle with respect to its rest position, *p* is its conjugate momentum, the *s* superscript stands for “standard,” and *g*(*t*) is the coupling function. Note that, in our approach, we keep *g*(*t*) general. By rescaling the variables as $q=m\Omega x$ and $p\u2192(1/m\Omega )p$, we rewrite the previous equation as

We shall define $H1s(t)=(1/m\Omega )g(t)a\u2020aq$ as the interaction part, and the free part as $H0=(\Omega /2)(p2+q2)+\u210f\Omega F(a\u2020a+(1/2))$.

We seek to express the resulting oscillator-state-independent phase accumulated by the light field at the end of the interaction. The fact that this phase term does not depend on the oscillator state is important as we seek to test modifications to the commutation relations. With no knowledge of the true commutation relations, the initial (ground) state of the mechanical oscillator is ill-defined. In the interaction (or Dirac) picture, operators *A* are transformed as $A\u0303(t)=eiH0t/\u210fAe\u2212iH0t/\u210f$. The unitary evolution of the system can be written as $\rho \u0303(t)=U\u0303(t)\rho \u0303(0)U\u0303\u2020(t)$, and since its time derivative must be consistent with the Liouville equation $(d/dt)\rho \u0303(t)=\u2212(i/\u210f)[H\u03031(t),\rho \u0303(t)]$, this requires the unitary propagator to satisfy $(d/dt)U\u0303(t)=(\u2212i/\u210f)H\u03031(t)U\u0303(t)$. Since the interaction Hamiltonians do not necessarily commute at different times, the propagator solving the equation takes the form $U\u0303t=exp\u2009(\Theta (t))$, where $\Theta (t)=\u2211k=1+\u221e\Theta k(t)$ is given by a Magnus series.^{26}

Because the first order unequal time commutators are scalars, only the two first terms in the Magnus expansion are non-vanishing. Explicitly, we have

With those two terms, we have exactly solved the dynamics, as the interaction Hamiltonian at any time commutes with the commutator, all following terms of the Magnus expansion are zero. To lighten the notations, let us rewrite the interaction term as $H\u03031s(t)=a\u2020a(g\u2032(t)q+\u2009g\u2033(t)p)/m\Omega $, where $g\u2032(t)=g(t)\u2009cos\u2009(\Omega t)$ and $g\u2033(t)=g(t)\u2009sin\u2009(\Omega t).$ The first term in the Magnus expansion then reads

where $G\u2032=\u222b0tg\u2032(t1)dt1,\u2009G\u2033=\u222b0tg\u2033(t1)dt1$, and the second term reads

where we have applied the standard commutation relation $[q,p]=i\u210f$. Hence, we find that the propagator takes the form

where $F(t)=12\u222b0tdt1(g\u2032(t1)G\u2033(t1)\u2212g\u2033(t1)G\u2032(t1))$ gives a phase term that is independent of the mechanical oscillator's state and is quadratic in the photon number.

## III. THE GUP MODIFICATION FOR A SINGLE PARTICLE

We now introduce the GUP modification to the commutation relations. While it is possible to analytically work out the dynamics of the system by imposing $[x,p]=i\u210f(1+\beta p2)$, we instead work with the usual commutation relations but redefine the momentum as $p\u2192p(1+(\beta /3)p2)$, as in Refs. 15 and 27, which yields the modified commutation relations. The Hamiltonian, with reduced variables, now takes the form

We will not carry out any perturbative treatment of the coupling term *g*(*t*), but we simply expand the Hamiltonian in orders of *β* and keep the order 1 correction, assuming that higher order terms are negligible.^{28} We furthermore neglect possible corrections to the optical Hamiltonian. Let us also note that the extra interaction term due to the GUP correction changes the ground state, which is why we will focus on a state-independent nonstandard signature.

We shall now define $H1(t)=(1/m\Omega )g(t)a\u2020aq+(1/3)m\Omega 2\beta p4$ as the interaction term. In the same manner as the standard phase was found, we work in the interaction picture where the propagator satisfies the differential equation $(d/dt)U\u0303(t)=\u2212(i/\u210f)H\u03031(t)U\u0303(t)$. In order to compute the Magnus terms, one needs to express the different commutators between $H\u03031$ at different times, and while the *p*^{4} term does give rise to many non-vanishing higher-order commutators, to first order in the GUP correction *β*, the Magnus expansion only has five non-vanishing terms. We have explicitly

where usual harmonic oscillator dynamics give $q\u0303(t)=cos\u2009(\Omega t)q+\u2009sin\u2009(\Omega t)p$ and $p\u0303(t)=cos\u2009(\Omega t)p\u2212sin\u2009(\Omega t)q$. Using standard commutation relations, one establishes the following commutator expressions:

To order 1 in *β*, it is clear that the first commutator $[H\u03031(t1),H\u03031(t2)]$ will take the form $C+\beta p\u03033+O(\beta 2)$, where *C* commutes with the interaction Hamiltonian. Therefore, the next order nested commutators are all at least order 1 in *β*. After some calculations, see Appendix A for details, one arrives at the following form of the fifth order commutator:

where $(t4\u2194t5)$ signifies that there is an extra term that is identical up to a swap between two time arguments.

The fifth order Magnus term, which is an integral of such fifth order nested commutators (up to some permutations of time indices), is particularly interesting as it is independent of the mechanical-oscillator operators $q\u0303(t)$ and $p\u0303(t)$. In order to further work out its form, we note that it is a quintuple nested integral of 22 nested commutators^{29} as shown explicitly in Appendix B. Simplifications are possible as 18 of those commutators are of the form $[[H\u03031(ti),H\u03031(tj)],[H\u03031(tk),[H\u03031(tl),H\u03031(tm)]]]$, and looking at Eqs. (A3) and (A4), we see that they are at least quadratic in *β*. To first order in *β* we are left with the four terms

where we have defined the nested integral operator $\u222b(5,t)dt5:=\u222b0tdt1\u222b0t1dt2\u222b0t2dt3\u222b0t3dt4\u222b0t4dt5$, the set *S* which is made up of the four permutations $(\sigma 1,\sigma 2,\sigma 3,\sigma 4)=((54321),(15423),(14325),(15324))$ and the associated coefficients $\lambda \sigma 1=\lambda \sigma 3=\lambda \sigma 4=\u2212(1/30)$ and $\lambda \sigma 2=(2/15)$. Combining the previous equation with Eq. (11), one arrives at the following expression that holds for arbitrary coupling functions *g*(*t*), to first order in *β*:

Similar to the pulsed regime considered in Ref. 22, this fifth Magnus term, which is proportionnal to the GUP correction *β*, contributes to the propagator as a highly non-linear phase factor depending on the fourth power of the photon number.

This is to compare with the standard case of unmodified commutation relations, where the analogous oscillator-state-independent phase factor $exp\u2009(\u2212(i/\u210f)m(a\u2020a)2F(t))$ in Eq. (7), is quadratic in the photon number. Our formula is valid for any optomechanical coupling function *g*(*t*). In Sec. IV, we shall analytically apply it to the case where the coupling function consists of four pulses, as was considered in Ref. 22.

## IV. PHASE TERM IN THE PULSED REGIME

We now make use of our general formula for a single particle, with a specific and relatively simple coupling function. To do so, it is useful to first define $f(t1,\u2026,t5)=\u220fj=14\u2009cos\u2009(\Omega (t5\u2212tj))$. In the context of pulsed interaction, the coupling term is well approximated by a sum of Dirac distributions $g(t)=\lambda \u2211i=03\delta (t\u2212\tau i)$. The Dirac delta interaction leaves the mechanical oscillator's state fixed over the duration of the pulse, which is identical to the approximation used for pulsed optomechanics (up to a global phase on the light). Here, a succession of four Dirac pulses at instants $(\tau i)0\u2264i\u22643$ is considered, where $\lambda \u22650$ is a coupling strength.

After some calculations shown in Appendix C, we find

where $H(a,b,c,d,e)$ is the succession of Heaviside functions $H(a\u2212b)H(b\u2212c)H(c\u2212d)H(d\u2212e)$.

In the particular case where the interaction instants are $\tau k=t0+k\theta =t0+k(\pi /2\Omega )$, where *θ* is a quarter of the oscillator period, and for times $t>t0+3\theta $ after which all four matter-field interactions have taken place, one finds the closed form result

All higher order Magnus terms vanish; hence, the solution obtained from the $\Theta i,1\u2264i\u22645$ is exact, to order 1 in *β*. We observe that this fifth-Magnus term is a constant of time for $t>t0+3\theta $, which is expected as the full system undergoes free evolution. Similar to what was found in Ref. 22, this mechanical-state-independent phase factor increases as the fourth power of the optical intensity $(a\u2020a)$ and the interaction strength *λ*. In Secs. V and VI, we will be studying how this nonstandard phase factor scales with the number of mechanical oscillators.

## V. GENERALIZATION TO THE MANY-BODY CASE

In view of the observations made in Refs. 22 and 24, the GUP correction *β* is expected to be linearly or quadratically suppressed by the size of the considered systems when they are uncorrelated or rigid. However, it may well be that in some intermediate regime, where a many-body system is neither rigid nor a collection of uncoupled particles, a nonstandard signature could scale superquadratically and thus overcome the suppression of *β*. We have taken our nonstandard signature to be the fifth order Magnus term to first order in *β*, which is an oscillator-independent phase factor. Here, we investigate the form of this nonstandard phase factor for many-body systems.

Let us extend our analysis to the case of *N* harmonically coupled identical trapped particles. The total Hamiltonian to first order in *β* can be written as $H(t)=H0+H1(t)$, where the free Hamiltonian is

However, our approach does not require a specific potential term and holds for general potentials of the form $(m/2)xThx$, where *h* is a symmetric matrix (Hessian of the potential). The interaction Hamiltonian is

where Ω_{c} is the coupling frequency between neighboring oscillators. In the same way as before, the time propagator in the interaction picture satisfies the differential equation $(d/dt)U\u0303(t)=\u2212(i/\u210f)H\u03031(t)U\u0303(t)$. Thus, the propagator is expressed as a Magnus series $U\u0303(t)=exp\u2009(\u2211k=1+\u221e\Theta k(t))$. We are lead to expressing nested self-commutators of the interaction Hamiltonian. Expressing the interaction Hamiltonian (17) in the interaction picture, one arrives at an expression similar to Eq. (12) for the fifth order Magnus term:

where to first order in *β* one has

In order to express the nested commutator in general, we claim that $[x\u0303i(t1),p\u0303j(t2)]:=Cij(t1\u2212t2)$ is a *c*-number. Indeed, as shown in details in Appendix D, $Cij(t1\u2212t2)=i\u210f\u2211k=1NOikO\u2032jk\u2009cos\u2009(\omega k(t2\u2212t1))$, where *O* and $O\u2032$ are some scaled normal mode transformation matrices defined by $O=(1/m)P\omega \u22121/2$ and $O\u2032=mP\omega 1/2$, where the orthogonal diagonalizing transformation *P* of *h* satisfies $diag((\omega i2)1\u2264i\u2264N)=PThP$. Note that *m* is a scalar and $\omega =diag((\omega i2)1\u2264i\u2264N)$ is a diagonal matrix.

One can now express the nested commutators as

and Eq. (19) can be cast as

where $Dj(t,t\u2032)=\u2211i=1Ngi(t)Cij(t\u2212t\u2032)$. Hence, the fifth Magnus term is found to read

where one can reduce the problem to the knowledge of the normal modes and frequencies by the explicit formula $Dj(t,t\u2032)=i\u210f\u2211i=1Ngi(t)\u2211k=1NOikO\u2032jk\u2009cos\u2009(\omega k(t\u2212t\u2032))$.

We have thus found an analytical form for the mechanical-state-independent phase factor. As in the single oscillator case, this factor depends on the fourth power of the optical intensity; however, the expression now involves the mechanical properties of the body described by its normal mode transformations and associated frequencies. To no surprise, if the potential terms correspond to coupled trapped particles, one recovers Eq. (13) by setting *N *=* *1, for which the normal mode transformations are trivial.

Finding how $\Theta 5(t)$ should scale with the number of particles *N* generally is not evident. In what follows, we will consider once again a pulsed regime where four pulses of light interact with the system, which will consist of many copies of a mechanical oscillator that are coupled to the neighboring oscillators, and show numerically that there can be superlinar scaling.

## VI. EXTENDING THE FOUR PULSE SCHEME

We seek to use our general formula given by Eq. (22) in the case where four pulses of light are sent through a lattice of *N* coupled oscillators (or sites). That is, we consider coupling functions of the form $gk(t)=\lambda \u2211i=03\delta (t\u2212(t0+iT+(k\u22121)\tau ))$, where *λ* is a coupling strength, *T* is the time separating two pulses for a given site, and *τ* the time taken by the pulse to travel from one site to the next one. We shall work under the assumption that $T\u2265N\tau $, which means that the pulse is re-injected only once the it has finished interacting with all the lattice sites. The setup and the interaction functions are illustrated in Fig. 1.

After calculations similar to those carried out in the context of a pulsed regime with a single particle, see Appendix E for details, one arrives at the following expression for the fifth order Magnus term:

where $\phi \nu 1\nu 2\nu 3\nu 4(t1,\u2026,t5)=\u220fs=14\u2009cos\u2009(\omega \nu s(ts\u2212t5))$ and $\u2200r\u2208{1,2,3,4},\u2009\theta r:=t0+\alpha rT+(ir\u22121)\tau $.

By fixing *N *=* *1 one naturally recovers Eq (14). In the case where the coupling between oscillators is zero, then the $O,O\u2032$ normal mode transformation matrices are diagonal, and one is left with the single oscillator phase term multiplied by *N*. Hence, in the uncoupled setting, one expects a linear scaling of the phase contribution coming from the fifth Magnus term with the number of oscillators, as confirmed by numerics.

As shown in Fig. 2, numerics reveal the presence of superquadratic scaling, at least up to *N *=* *5, for coupling frequencies Ω_{c} equal to a tenth or a half of the trap frequency that we chose to be $\Omega =2\pi \xd7105\u2009Hz$ to facilitate comparison with Ref. 22. The plotted quantity is an absolute phase factor $|\Phi (N)|$ for $N\u2208{1,2,3,4,5}$ defined by $\Theta 5(t)=(\u2212i/\u210f)(4!/3)(\beta /m)(a\u2020a)4(\lambda 4/30)\Phi (N)$. In general this is a complicated factor, depending on the coupling *g*(*t*), the trapping frequencies of each lattice sites, and the coupling frequencies between them. We have assumed the form given by the last equation, i.e., we work in pulsed regime, and the pulses are separated by a quarter trapping period which is assumed to be the same for all sites, and we have further assumed $\tau =T/(2N)$.

For coupling frequencies comparable or higher to the trap frequency, we see in Fig. 3 that there remains an advantage; however, for the considered coupling frequencies, the absolute phase factor is not an increasing function of the number of sites.

## VII. CONCLUSION

We have revisited modern approaches to indirectly test non-classical features of spacetime, namely, an optomechanical protocol that aims to detect highly non-linear phase terms that should arise as a result of modifications of commutation relations. Although rigid macroscopic systems are not viable probes for Planck scale effects, as the GUP corrections to the commutation relations of center of mass variables of such systems are expected to scale as the squared inverse of the number of constituents, the case of correlated many-body systems was left unexplored.

In this paper, we have considered such systems in general and, to first order in GUP corrections, we have derived the full dynamics for completely general light–matter coupling functions, and found explicit forms for the fifth order Magnus term that consists of a phase factor that is independent of the mechanical oscillator's state. Many parameters of correlated many-body systems could be tuned, such as the retardation time between two lattice sites, or the spatial width of the optical pulse. Limiting ourselves to a pulsed regime, we have shown that superquadratic scaling is already possible up to *N *=* *5 subsystems by tuning only the coupling frequency between lattice sites.

This provides new perspectives for controlled many-body quantum dynamics and may help to re-affirm finely tuned macroscopic probes as potential amplifiers for testing quantum gravitational effects. In addition to showing how macroscopic probes may be relevant to test the Planck scale physics, the general calculations laid out in this work may benefit further investigations and tests of fundamental physics using optomechanical systems.

## ACKNOWLEDGMENTS

We thank D. Jennings for helpful comments. This work is supported by the Korea Institute of Science and Technology (KIST) Open Research Program. H.K. is supported by the Korea Institute for Advanced Study (KIAS) Individual Grant No. CG085301 at Korea Institute for Advanced Study. I.P. was supported by The Branco Weiss Fellowship—Society in Science and the European Research Council (ERC) under Grant No. 742104. This work was funded by the Engineering and Physical Sciences Research Council (ESPRC) Centre for Doctoral Training in Controlled Quantum Dynamics, and QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union's Horizon 2020 Programme and KIAS visiting professorship.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: DETAILED EXPRESSIONS FOR THE NESTED COMMUTATORS

Recall that the interaction Hamiltonian in the interaction picture reads

Then, to first order in *β* one has

Using the commutator expressions (10), the previous equation takes the explicit form

Since the first term commutes with $H\u03031$, the next nested commutator to first order in *β* takes the simple form

The following term is then found to be

from which one finds Eq. (11).

### APPENDIX B: THE GENERAL FIFTH ORDER MAGNUS TERM

Using the shorthand $H\u0303i:=H\u03031(ti)$, it was shown explicitly^{29} that

As stated in the main text, with the interaction Hamiltonian given by Eq. (9), only the four first terms are linear in the GUP correction *β*. Indeed, the first order commutators $[H\u0303i,H\u0303j]$ are of the form $C+\beta p\u03033$, where *C* always commutes with the interaction Hamiltonian. The second order commutators given by Eq. (A4) are at least linear in *β*. Hence, all 18 commutators of the form $[[H\u0303i,H\u0303j],[H\u0303k,[H\u0303l,H\u0303m]]]$ are at least quadratic in *β*.

### APPENDIX C: DETAILED DERIVATION OF THE PHASE TERM IN THE PULSED REGIME

We begin with Eq. (13)

which we rewrite as

The *λ _{k}* are known. We seek to express the

*T*terms when the interaction

_{k}*g*takes the form $g(t)=\lambda \u2211i=03\delta (t\u2212\tau i)$ with no further specification for

*τ*. To this end, we introduce $f(t1,\u2026,t5):=\u220fj=14\u2009cos\u2009(\Omega (t5\u2212tj))$. Then

_{i}Showing explicit calculations for the first term $\sigma 1=(54\u2009321)$, we have

Integrating the *δ* distributions yields

In a similar fashion, from $\sigma 2=(15\u2009423)$, we obtain

For the third permutation $\sigma 3=(14\u2009325)$, one finds

Finally, from $\sigma 4=(15\u2009324)$, it follows that

Summing up these expressions, noting that $\lambda 1=\lambda 3=\lambda 4=\u2212(1/30)$ and $\lambda 2=(4/30)$, one arrives at

Let us now specify $\u2200k\u2208{0,1,2,3},\u2009\tau k=t0+k\theta $, where $\theta =(\pi /2\Omega )>0$ and assume $t>t0+3\theta $ then we find

Since $f(t0+i1\theta ,t0+i2\theta ,t0+i3\theta ,t0+i4\theta ,s)=\u220fk=14\u2009cos\u2009(\Omega (s\u2212t0)+ik\pi /2),$ the affine change of variable $u=\Omega (s\u2212t0)$ yields

### APPENDIX D: DERIVATION OF THE $[x\u0303,p\u0303]$ UNEQUAL TIME COMMUTATOR

We prove our claim about $[x\u0303i(t1),p\u0303j(t2)]$ being a *c*-number and derive the explicit form of $Dj(t,t\u2032)$. To this end, one can write the free Hamiltonian *H*_{0} as

where *h* is a symmetric matrix characterizing the potential term. One can define normal variables $(\varphi ,\pi )=(\omega 12PTmx,\omega \u221212PT(1/m)p)$, where $\omega 2=PThP=diag(\omega i2)$ is the diagonal form of the *h* matrix, consisting of the squares of the normal frequencies, and where *P* is the orthogonal transformation matrix. Let us stress that *m* is a scalar and *ω* is a diagonal matrix. With those normal variables, the free Hamiltonian takes the simple form

We furthermore note that the transformation is canonical as $[\varphi i,\pi j]=i\u210f\omega i/\omega j\u2211k,l=1NPkiPlj\delta k,l=i\u210f\delta ij$ by the fact that $PTP=IN$. For simplicity, we will rewrite the variable transformations as $x=O\varphi $ and $p=O\u2032\pi $, where *O* and $O\u2032$ are orthogonal matrices. Using the shorthand notation $U0(t):=e\u2212iH0t/\u210f$, one can now establish that

We have $\varphi \u0303k(t1)=cos\u2009(\omega kt1)\varphi k+sin\u2009(\omega kt1)\pi k$ and $\pi \u0303\u2113(t2)=cos\u2009(\omega \u2113t2)\pi \u2113\u2212sin\u2009(\omega \u2113t2)\varphi \u2113$, from which it follows that

which proves that $Cij(t1\u2212t2)$ is a *c*-number.

### APPENDIX E: EXTENDING THE OPTOMECHANICAL PROTOCOL SCHEME

One can make Eq. (22) more explicit by writing out

where $fi1i2i3i4j(t\sigma 1,\u2026,t\sigma 5)=\u220fs=14Cisj(t\sigma s\u2212t\sigma 5)$. We have $Cisj(t\sigma s\u2212t\sigma 5)=i\u210f\u2211\nu =1NOis\nu O\u2032j\nu \u2009cos\u2009(\omega \nu (t\sigma s\u2212t\sigma 5))$. Thus, one can write

where $\phi \nu 1\nu 2\nu 3\nu 4(t\sigma 1,\u2026,t\sigma 5)=\u220fs=14\u2009cos\u2009(\omega \nu s(t\sigma s\u2212t\sigma 5))$. From this, one can express the fifth Magnus term as

We may write

where $\theta (\alpha r,ir):=t0+\alpha rT+(ir\u22121)\tau $ for convenience. We shall lighten the notation to *θ _{r}* bearing in mind that it is a function of

*α*and

_{r}*i*. The task is now reduced to expressing the integral

_{r}After similar calculations as carried out in Appendix C for the single oscillator driven by pulsed light, one arrives at

where the coefficients associated with the permutations are $(\lambda 1,\lambda 2,\lambda 3,\lambda 4)=((\u22121/30),(4/30),(\u22121/30),(\u22121/30))$ and the integration domains are given by oriented intervals (and unions of)

Thus, one can rewrite the fifth Magnus term as

Note that with no coupling between the oscillators, $Ois\nu sO\u2032j\nu s=\delta is,\nu s\delta j,\nu s$ so we are left with

The vector of normal frequencies is simply $\omega T=\Omega (1,1\u2026,1)$ and $\theta k=t0+\alpha kT+(j\u22121)\tau $, so the Heaviside functions only depend on the *α* coefficients.

The functions $\phi jjjj$ are thus identical and the sum over *j* simply gives a multiplication by the number of lattice sites *N*. That is, one arrives at *N* times the Magnus term given by Eq. (C9). The absence of coupling between the oscillators is equivalent to repeating the single oscillator experiment a number of times, from which one does expect a linear increase in the nonstandard signature picked-up by the light.

## References

*β*

_{0}can be interpreted as effectively reducing the value of $\u210f$ up to some momentum cutoff.

*β*come with higher powers of the operators. The assumption that higher order corrections are negligible assume some reasonably localized state in the (

*q*,

*p*) phase space, which is the case for low energy states.