Enhancing a slow and weak optomechanical nonlinearity with delayed quantum feedback

A central goal of quantum optics is to generate large interactions between single photons so that one photon can strongly modify the state of another one. In cavity optomechanics, photons interact with the motional degrees of freedom of an optical resonator, for example, by imparting radiation pressure forces on a movable mirror or sensing minute fluctuations in the position of the mirror. Here, we show that the optical nonlinearity arising from these effects, typically too small to operate on single photons, can be sufficiently enhanced with feedback to generate large interactions between single photons. We propose a protocol that allows photons propagating in a waveguide to interact with each other through multiple bounces off an optomechanical system. The protocol is analysed by evolving the full many-body quantum state of the waveguide-coupled system, illustrating that large photon–photon interactions mediated by mechanical motion may be within experimental reach.

T he interaction between light and motion in cavityoptomechanical systems has enabled sensitive measurements of force and displacement, as well as quantum-optical control over the state of mechanical resonators 1 . Conversely, an optically coupled mechanical degree of freedom can lead to interactions between photons mediated by motion. Nonlinear optical effects such as wavelength conversion [2][3][4][5] , generation of squeezed light [6][7][8] and electromagnetically induced transparency 9,10 are manifestations of this optical nonlinearity and have been demonstrated in recent experiments. Nonetheless, a central goal of quantum optics and information is to generate nonlinearities that are large at the single-photon level 11 . It has been shown theoretically that it is possible to generate singlephoton effective optical nonlinearities with optomechanics, but only with system parameters well outside of our current experimental reach [12][13][14][15][16][17][18] . The essential difficulty is that for a mechanically induced optical nonlinearity to act on the incident light, the response of the mechanical system must be fast compared to the amount of time the photon spends inside the optical cavity. Moreover, the mechanical system should display a large response to the force induced by a single photon-so that the cavity properties change significantly due to the presence of the photon-requiring a highly compliant mechanical system. These two contradicting needs, fast response and large compliance, are at odds with each other and make strong effective photon-photon interactions extremely difficult to achieve in realistic cavity-optomechanical systems.
In addition to the difficulty of generating the large Kerr nonlinearities needed, it has been pointed out that in principle, subtle effects due to the multimode nature of a propagating light field make implementation of high-fidelity gates operating on flying photons using large instantaneous Kerr nonlinearities impossible 19 .
Recently, it has been proposed that optomechanical nonlinearities can be enhanced in a fast-cavity system by parametric amplification of mechanical motion 16 leading to a quantum gates between photons of different frequency. We propose a different way to make strong photonic interactions from weak optical nonlinearities. In this article, we show that by introducing a coherent delay to a cavity-optomechanical system with a small optomechanical coupling and slow response, large photon-photon interactions between two temporally separated photon modes propagating in an optical fibre can be achieved. Time-delayed coherent feedback has been discussed for Gaussian states in other contexts 20,21 . Recently a new method based on matrix product states (MPS) 22 has been developed that greatly facilitates a full quantum analysis of such interactions. Here, we demonstrate using both a semiclassical argument and full quantum simulations taking into account the propagating quantum field with its many degrees of freedom, that a quantum phase gate between two temporal modes of the field can be implemented with a cavityoptomechanical system that is in the bad-cavity regime. The coherent delay allows a slow optomechanical system to induce a large effective interaction between the temporally separated photons, greatly relaxing the optomechanical system requirements needed to implement such a gate. The nonlocal nature of the gate also sidesteps a key assumption in the aforementioned impossibility arguments 19,23,24 .

Results
System and phase gate protocol. We consider an optomechanical system coupled to a long waveguide with an end mirror, so that photons can propagate back and forth inside the waveguide, and interact repeatedly with the mechanical resonator by entering the optical cavity. The cavity-optomechanical system is composed of an optical resonator at frequency o o with annihilation operatorâ , coupled to a mechanical resonator with frequency o m and annhililation operatorb. The initial state of the waveguide consists of two temporal modes of the light field with Gaussian profile, each with extent t. We assume that 1=t ( k, where k is the optical loss rate, so that each photon can fully enter the cavity and exert a significant impulse onto the mechanical system. During this process, the photon also obtains an uniform phase shift caused by the internal state of the mechanical system. Each impulse is assumed to be nearly instantaneous on time scales relevant to the mechanical oscillator's motion, that is, o m ( 1=t. Taken together, these conditions imply that we are operating in the bad-cavity regime o m ( k. The slow internal dynamics of the mechanical system as well as the coherent time delay lead to an effective nonlinear photonphoton interaction that is nonlocal in time and also lead to a build up of entanglement between the two temporal modes inside the waveguide. Remarkably, we find that a quantum controlled phase gate (CPHASE), where |00i, |01i and |10i are left unchanged and |11i-À |11i can be implemented in this system.
The protocol for the CPHASE gate between two temporal modes of the light field is shown in Fig. 1. The centre frequencies of these two modes are chosen to be at the cavity frequency and each have a bandwidth much smaller than the cavity linewidth. They are also separated by a time T m =4 ) t, where T m ¼ 2p/o m is the period of the mechanical oscillation. Our system is similar to and inspired by pulsed optomechanical experiments 25,26 with the distinction that no measurement occurs in our feedback networkthe photons are fed back to the system coherently. The choice of delay time is essential for disentangling the mechanical system from the photons at the end of the protocol. The mechanical system must return to its initial state, the ground state, regardless of the state of the input temporal modes, since any residual entanglement between the temporal modes and the mechanical system will reduce the fidelity of the CPHASE gate. As shown in Fig. 1 and is justified below, this condition is satisfied for the temporal mode spacing we have chosen.
Semiclassical model. Here we use a simplified semiclassical model to understand some of the behaviours of this system. The optomechanical system's Hamiltonian is given bŷ The radiation pressure force on the optomechanical system ðF RP ¼ À ' g 0â wâ =x zp Þ integrated over the interaction time of the photon with the cavity t ( o À 1 m À Á causes a rapid change in the momentum of the mechanical system by Dp RP ¼ À ' g 0 R t dtâ y ðtÞâðtÞ=x zp . An input fieldĉðtÞ from the waveguide is incident on the cavity at time t. If we assume that the cavity can be eliminated adiabatically (k is large), the field in the cavity and waveguide may be related asâðtÞ % 2ĉðtÞ= ffiffiffi k p , so after the interaction with the propagating photon, an impulse of Dp RP ¼ À 4' g 0 R t dtĉ w ðtÞĉðtÞ=x zp k is imparted onto the mechanical system. We define a photon number operator n k ¼ R k dtĉ w ðtÞĉðtÞ, which counts the number of excitations in the kth temporal mode. The interaction between the kth temporal mode and the mechanical system is given by the unitary operator with positionx¼x zp ðb w þbÞ. The mechanical free evolution operator isÛ t ¼e À io m b w bt . Our protocol can then be compactly stated asÛ protocol C j i wg 0 j i m , witĥ This sequence of interactions is shown schematically in Fig. 1. We can now calculate the result of these operations on the joint state of the optomechanical system and photonic waveguide. We take state of the photonic waveguide to be initialized as |Ci wg ¼ |ji 1 |ki 2 , the state with j photons in the first temporal mode and k photons in the second temporal mode. The mechanical mode and optical cavity are assumed to be in their ground states. Then where f jk ¼ 0 for jk ¼ 00, 01, 10 and f jk ¼ f 1 for jk ¼ 11 with According to this simplified model, obtaining a total phase shift of p can be accomplished either by going to a large enough coupling g 0 /k, or by running the protocol multiple times, N rep Ep/f 1 , for a smaller g 0 /k so that a total phase shift builds up over several bounces. This requires controlling the number of bounces, which can be accomplished, for example, by using high extinction Mach-Zehnders that are being developed for photonic quantum information processing. To distinguish between these two approaches, it is important to balance the losses incurred over multiple bounces with the deleterious effect of using a larger coupling and fewer bounces. For losses, we consider that in any real system, in addition to the coupling between the cavity and the waveguide, there are other intrinsic optical cavity losses characterized by a loss rate k in leading to a probability k in /k that a photon is lost on each bounce. According to the simplified semiclassical model, since the chance of photon loss in multiple runs is proportional to N rep k in /k, and N rep pk 2 , it is always advantageous to make k as small as possible. This conclusion however neglects effects that arise due to strong coupling. One of these effects is the large optical frequency shift caused by the mechanical displacement induced by the first photon, which prevents the second photon from fully interacting with the system. This prevents a perfect erasure of the information about the photons from the mechanical oscillator and causes the residual photon-phonon entanglement. A model described in the Methods section shows that interaction of the waveguide state |1i 1 |1i 2 with the optomechanical system leads to a final state e if |1i 1 |1i 2 |b r i m after one repetition of the protocol, where the mechanical system is in a coherent state |b r i and b r p(g 0 /k) 5 . Other input states do not cause a change in the mechanical state. This residual phonon occupancy for the |11i input state leads to a reduction in gate fidelity on the order of 0jb r h i j j 2 ¼e À b r j j 2 which can be made very small by going to smaller coupling and a larger number of bounces. Other effects caused by the finite extent of t should also be considered carefully. For example, a photon wavepacket can obtain a non-uniform position-dependent phase shift due to frequencydependent phase response of the cavity, causing the state of the field to no longer be in the same temporal mode. Such an effect is analysed in more details in the Method section and can be made smaller by making kt larger. Also a mechanical system can have a finite position shift during a bounce of DxEDp RP t/m, which causes another non-uniform phase shift on the order of 4g 0 Dx/kx zp ¼ f 1 o m t for a single-photon input state in a given temporal mode. This effect can be reduced by making o m t smaller. Both of these types of imperfection cause the state of the electromagnetic field to move out of the subspace spanned by |ji 1 |ki 2 |0i m and are difficult to capture quantitatively in an analytical form. Estimates of fidelity given these effects are provided in the Methods section. However, quantitative calculations are important for understanding the realizability of the protocol given state-of-the-art experimental capabilities and therefore we turn to full quantum simulations to incorporate all these effects.
MPS simulation for quantum dynamics. The feedback network we are considering requires a long delay due to the slow dynamics of the mechanical oscillator. Recently methods for understanding dynamics of systems in such quantum feedback networks have been proposed 22,[27][28][29] . In addition, we are interested in understanding the evolution of the state of the photons in the waveguide and interactions induced between them by the optomechanical system. One approach to solving the full dynamics of the waveguide and system together is to discretize the waveguide into time steps Dt that are much smaller than any of the relevant dynamics of the system and numerically evolve what is now effectively the interaction between a one-dimensional (1D) chain of harmonic oscillators and the system 22 . In principle, keeping track of the state of such a 1D chain is daunting due to the exponentially large Hilbert space that scales as O d Tm=2Dt À Á , where d is the truncated dimension of each time-bin's Fock space. Luckily, the states of the waveguide we are considering have far less entanglement than general states in the full Hilbert space, so  Figure 1 | Protocol for the CPHASE gate. Two temporal photonic modes in the waveguide are separated by a time interval of T m /4. The modes can be in states |0i and |1i. First, temporal mode 1 interacts and is entangled with the mechanical oscillator as |0i and |1i will lead to different changes in the mechanical momentum and cause a state-dependent change in the mechanical state. After the mechanical system evolves for T m /4, the second photon mode interacts with the mechanical oscillator and the mechanical system's state become entangled with both photons. After the mechanical system evolves again for another T m /4, the first photon temporal mode interacts again with the mechanical system causing its state to be disentangled from the mechanical oscillator. After another T m /4 of evolution, the second photon mode comes back again and the mechanical system is decoupled from both temporal modes since it goes back to its ground state regardless of the initial states of the photon modes. The blue lines signify the entanglement in the system and on the right the statedependent evolution of mechanical state is shown in phase space. efficient methods for storing and evolving the states can be utilized 22,[30][31][32][33][34][35] . Our full quantum simulation is based on the MPS representation of the quantum field. MPS as well as the related time-dependent density matrix renormalization group techniques have already been well-established in condensed matter physics for simulating 1D quantum many-body systems.
To implement the MPS method for the continuous quantum field of the waveguide, we discretize time in small steps Dt and define operatorĉ n ¼ 1 ffiffiffi ffi Dt p R t n t n À Dtĉ ðtÞdt for the nth time-bin, where t n ¼ nDt, n is an integer andĉðtÞ is the field operator in time domain. Using the commutation relation ½ĉðtÞ;ĉ w t 0 ð Þ¼d t À t 0 ð Þ, it is straightforward to verify that ½ĉ n ;ĉ w m ¼d nm , which means the time-bins can be interpreted as independent harmonic oscillators. The state of this quantum many-body system can be represented in the canonical MPS form ref. 30. A spatially distributed single-photon state of the waveguide is then 1 We consider an initial state where the optomechanical system's optical and mechanical modes are both in their ground states, and consider two temporal modes of the photon with annihilation operatorsÂ 1 andÂ 2 separated by a time interval T m /4. We identify the states jk j i m for j, k ¼ 0 or 1, and assume an initial state of the whole system Note that the states |jki are to an extremely good approximation orthogonal given a temporal separation that is much larger than their width. As shown in the Methods section, initially the optomechanical system is at the first site of the MPS. To evolve the many-body state, we sequentially update the MPS by applying the unitaryÛ n ¼ exp À iĤ S Dt þ ffiffiffiffiffiffiffiffi kDt pâĉ w n Àâ wĉ n À Á À Á as a local gate on the nth time-bin and the optomechanical system. Then a swap gate is used to permute the order of them so that U n þ 1 can be applied in the same way as U n . The swap gate is used because it is more convenient to update the MPS in a local way by applying the gates only to the nearest neighbour sites-it does not represent a physical evolution of the many-body state. This process effectively simulates a discrete representation of the quantum stochastic schrödinger equation 22 and is continued until the optomechanical system reaches the last site of the MPS.
In contrast to ref. 22 where the number of sites in the MPS is proportional to the total simulation time, here our waveguide is modelled as a finite number of time bins corresponding to the feedback waveguide length. After every T m /2 time interval, corresponding to the total optical path length for a round trip, the optomechanical system has interacted with each time bin and reaches the last site of the MPS. In the next time step, due to the reflection off the far mirror, the system interacts with the first time bin again. This is accomplished by moving the system back to the first site of MPS via a series of swap gates. At this point, the evolution can be continued as before to simulate interaction after one bounce. This process is repeated 2N rep times for N rep runs of the protocol.
Entanglement entropy and phase calculation. MPS provides a convenient way to extract information about entanglement entropy of a quantum system. In the MPS representation we decompose the state vector of the N-site system, c i 1 ... i N , into a product of tensors G ½ki N a k À 1 a k and vectors l ½k a k for k ¼ 1yN, such that a bipartition of the state at bond l can be written as a Schmidt . This allows us to calculate the entanglement entropy between the two halves of the system by simply reading off the values of the l vector at bond l, and calculating S l ¼ À P a l l ½l2 a l log 2 l ½l2 a l . The black line in Fig. 2a is a plot of S l as a function of lDt/T m for the initial state |Ci i , which is essentially the bipartition entanglement entropy throughout the whole waveguide. The two peaks at the positions of the two photons correspond to the finite width photon excitations in the waveguide, since one excitation is spread across multiple local sites and detecting the photon at a particular bin tells us that there are no photons in the nearby bins. Stated more concisely, |1i f cannot be written as a product state of time-bin localized photon excitations. In addition to the localized peaks, two other interesting regions are that between the photons, and the region between the photons and the optomechanical system. We call the entanglement entropy between the photons and mechanical resonator S om and the entanglement entropy between the two photons S oo . When calculating S oo we are in fact finding the entanglement between one temporal mode of the waveguide and the rest of the system, including the other temporal mode as well as the optomechanical system. However, when S om is close to 0, that is, mechanical system disentangles from photons, S oo simply becomes the entanglement between two photon temporal modes. As expected, both S om and S oo are zero since the initial state can be written as a product state of separated temporal modes of the optical field in the waveguide and the optomechanical system, c.f. equation (6).
The entanglement entropy throughout the waveguide is plotted after a few subsequent runs of the protocol and shown in dashed lines in Fig. 2a. For the system parameters in Fig. 2a, S oo increases on each run of the protocol while S om remains close to zero. The evolution of S oo for more runs of the protocol is plotted in Fig. 2b and shows clear oscillatory behaviour. This can be attributed to the evolving phase shift of the |11i state with respect to the other states, since for a state c j i¼ 1 2 00 To verify that this is in fact the correct interpretation, we calculate the phase shift directly from the wavefunctions by taking the overlap between a vector 11 j i¼Â w  and the result of our simulation with non-zero g 0 after each run of the protocol. The reason for doing this instead of comparing with the initial unevolved state is that the temporal mode of the photon indeed changes but the distortion caused by a cavity can in principle be reversed using other linear optical components 36 . The total phase shift f increases linearly as N rep f 1 in Fig. 2d, which explains the oscillatory behaviour of entanglement entropy (in fact these phase shift gives exactly the same entanglement entropy as Fig. 2b). We compare the phase shift after one run of the protocol to the semiclassically predicted results f 1 ¼ 32(g 0 /k) 2 , and find good agreement at smaller g 0 /k as shown in Fig. 2c. The agreement becomes progressively worse at higher coupling. Finally, we calculate in the same way the phase shift for |00i, |10i and |01i components and find them all 0 for the values of g 0 /k in Fig. 2.
Fidelity of the CPHASE gate. In addition to calculating the phase, we verify that the state of the electromagnetic field after interaction with the optomechanical system remains within the subspace of states |jki, modulo the linear optical response of cavity. We calculate the fidelity F ¼ |h11|Ci| 2 after one and two bounces of the photons from the cavity. Nominally, after two bounces, that is, N rep ¼ 1, we expect F ¼ 1/4 for the input state in equation (6). In Fig. 3a a plot of the infidelity 1/4 À F against the interaction rate is shown for one and two bounces of the waveguide photons from the cavity. After one bounce, the significant entanglement between the mechanical system and waveguide photons leads to higher infidelity. After two bounces, that is, a full run of the protocol, the entanglement between the mechanics and waveguide photons is largely erased, causing an increase in the overlap F. Larger interactions cause a breakdown of this picture and lead to residual entanglement with the mechanical system and lower fidelity. For g 0 /k ¼ 0.1, we simulate a fidelity of 4FE0.9996 after two bounces, and 4FE99% after the N rep ¼ 10.
These simulations show that the photon wavefunction is not significantly distorted in the interaction, as would be expected for large instantaneous Kerr nonlinearities 19 .
To obtain a p phase gate, multiple bounces are required. In Fig. 3b, we show the simulated fidelity, F p ¼ |h11|Ci| 2 , and the required g 0 /k for phase gates implemented with different numbers of N rep . For phase gates with larger N rep , smaller g 0 /k are required; the gate fidelity is found to increase monotonically as the system becomes better described by the semiclassical model. Photon loss causes a reduction in the fidelity by a factor Z N rep , where Z is the probability of a photon being lost on a single repetition of the protocol. This finite photon loss causes the gate fidelity to be maximized at a finite N rep , with smaller N rep being favoured for higher losses Z.
Another signature of failure of the CPHASE gate is the presence of residual entanglement between the photons and the mechanical system after running the protocol, that is, a non-zero value for S om . In Fig. 3c,d, the evolution of S om and gate fidelity are plotted against g 0 /k and N rep . It is clear that increasing g 0 /k causes an increase in the residual entanglement and reduction in gate fidelity. The white points in Fig. 3d outline the relationship between N rep and g 0 /k needed to obtain a phase shift of p according to the quantum simulations. From these two figures, it is clear that the reduction in fidelity is largely due to residual photon-phonon entanglement, which can be approximated semiclassically. In the Methods section, we outline how the semiclassical model predicts a residual phonon occupancy giving qualitatively similar results to the quantum case, though the full quantum calculations are more forgiving in terms of obtainable phase shifts.
By studying the simulated fidelity for different parameters, we can make the phase gate conditions more precise. For N rep ¼ 4, we find that the bounds o m to0.3 and kt4200 are sufficient to prevent excess loss in fidelity. In our simulations the temporal width of the phonons t was chosen to be t ¼ 0.15/o m ¼ 1000/k.

Discussion
We have shown that strong photon-photon interactions can be obtained with an optomechanical system that is outside the strong coupling regime g 0 k=o 2 m 41 À Á . We estimate the effect of five sources of decoherence; (1) the mechanical coupling to the thermal bath, (2) the effect of other weakly-coupled mechanical modes and sources of phase fluctuations, (3) the intrinsic optical loss k in , (4) insertion loss and (5) the propagation losses in the long waveguide. The important parameter in understanding the mechanical decoherence (1) is the thermalization rate G m ¼ n b o m /Q m , where o m /Q m is the mechanical linewidth and n b is the thermal occupation at the mechanical frequency. In the high temperature limit, n b ¼ kT/:o m and G m ¼ kT/:Q m . As long as the G m ÂN rep T m ( 1, the chance of a phonon entering the system from the bath during the time that protocol is being run for remains small. This is equivalent to having f m Q m ) N rep kT=' . Such f m Q m products have been obtained at cryogenic temperatures 37 , and more recently at room temperature with silicon nitride membranes [38][39][40] . The effect of competing mechanical modes (2) has long been an issue in pulsed optomechanics experiments 25 . Recent experiments 26 have been successful at reducing the coupling to these modes to a few per cent of the primary mode by careful engineering of structures and positioning of the optical beam. In our case weak coupling to parasitic mechanical modes leads to dephasing and is considered more carefully in the Methods section. Sources (3)-(5) affect the protocol in a similar way, in the sense that they introduce a finite probability Z that a photon can be lost to the environment during the execution of the gate. We assume throughout this discussion that total loss rate k is composed of the intrinsic and extrinsic parts k ex þ k in , with k ex ) k in . For the intrinsic losses in the cavity (3), this means requiring Z 2Nrep i % 1 À k in =k ð Þ 2Nrep % 1 À 4N rep k in =k be close to 1. This linear reduction in fidelity is seen in Fig. 3b and makes smaller N rep more optimal. Insertion losses (4) can be considered in an identical way. Finally (5), the chance of a photon being absorbed in the long delay is related to the attenuation length of the fibre (on the order of 0.15 dB km À 1 in the telecom C-band) and the required propagation path needed in the fibre, N rep T m v fibre . To achieve a more than 90% chance for the photon to survive propagation in the delay for the N rep ¼ 4 case, we require a mechanical frequency 4260 kHz. Together these parameters mean that N rep ¼ 4 and (g 0 , o m , k)/2p ¼ (31.3 MHz, 260 kHz, 173 MHz) are needed to implement the CPHASE gate with FE80%. Currently, these parameters, though significantly easier to obtain than the large g 0 /k and good cavity limit koo m , remain outside of our experimental reach. Nonetheless, progress in experimental quantum optomechanics, and techniques for enhancing coupling by stacking low frequency membranes 17 , are expected to bring us to the required parameters in the coming years.
In conclusion, we have demonstrated that a CPHASE gate between propagating photons can be implemented with the slow and weak nonlinearity of an optomechanical resonator under the correct quantum feedback conditions. In addition to opening up a part of the optomechanical parameter space to quantum experiments, the approach shows that quantum feedback has the potential to enhance quantum nonlinear phenomena. In the future, it is interesting to consider the evolution of entanglement in quantum feedback networks using the MPS method, to better understand how interesting quantum many-body states with nontrivial correlations can be generated.

Methods
Here we will go through some details of the matrix product state method, giving an explicit algorithm to decompose a single-photon state with wide temporal extent into MPS form and also showing the equivalence between our simulation scheme and the dynamics of a closed waveguide system. We will also present results of fidelity simulations not included in the main text as well as a more detailed semiclassical analysis and the effect of parasitic mechanical modes.
Time-bin representation. The system-bath HamiltonianĤ¼Ĥ S þĤ B þĤ int , wherê H S ¼o mb yb þ g 0â yâby þb where we have definedĉ To calculate the evolution of the whole system, it is convenient to discretize time in small steps of Dt. Define the quantum noise increments in the time domain for the input field as:ĉ where t n ¼ nDt and n is an integer. It is straightforward to verify thatĉ n ;ĉ w m Â Ã ¼d nm from the commutation relationĉðoÞ; Þof the field mode, which means all the time-bins can be interpreted as independent harmonic oscillators. Thus the Hilbert space for the whole system is a tensor product space, including contributions from optomechanical system and quantum field in the waveguide, which can be modelled as a series of harmonic oscillators.
In the time-bin representation, the evolution of the system state in the nth time step is Initial state preparation. As usual, at t ¼ 0 the state of the optomechanical system and the quantum field in the waveguide are assumed to be completely uncorrelated and separable. The state of the waveguide is composed of two temporal modes that are spaced apart from each other. A large separation in time means that the state of the quantum field is very nearly the product state of two single-photon states in two independent temporal modes. A single photon in a temporal mode has a wavefunction given by 1 This can be rewritten in the time-bin representation as X n f nĉ More generally, the joint state of the waveguide with two temporal modes each in their 0 or 1 state can be described by the linear combination of state vectors Specifically, in our simulation the initial state is choosen to be C j i i ¼ 1 2 00 j iþ 01 j iþ 10 j iþ 11 j i ð Þ , as shown in the main text.
Matrix product state. For a general many-body state the following canonical decomposition always exists 30 : where the state is represented using a series of tensors G and vectors l. An important property of this canonical form is that for any k, l [k] is the Schmidt vector for a certain bipartition of the whole system into subsystem 1-k and subsystem k þ 1-n. For a quantum state with restricted amount of entanglement, the Schmidt vectors can be truncated by some threshold (10 À 4 -10 À 5 in our case-verifying that our results are not dependent on this threshold), which enables efficient classical simulation.
Single-photon state decomposition. Generally speaking, two things are required for the MPS simulation: preparing the many-body state |C(t)i in MPS form and updating this state step by step to evolve |C(t)i. In this and the next section, we consider these two aspects.
To decompose the initial state in an MPS form, a sequence of singular value decompositions are required. Note that for separated temporal modes each being in the 0 or 1 photon states, the dimension of Hilbert space for any of the time-bins can be chosen to be 2. In addition, we are operating in the polynomial O(N 2 ) subspace of the O(d N ) dimensional Hilbert space of the waveguide-only states |11? 0i, |101? 0i, ?|011? 0i, ?, |0? 11i are used since we ignore the possibility of two photons localizing at the same time-bin. Considering M to be the dimension of the truncated Hilbert space for the optomechanical system, the state of the full system would require O(N 2 M) parameters. The MPS representation takes advantage of this reduction in problem size in an implicit way. In addition, we note that though we are limited to this subspace due to the photon number conserving symmetry of the Hamiltonian and our selection of the initial state, our implementation of the MPS method is capable of simulating the dynamics for a more complex input state in an efficient and automatic way, since the Hilbert space is never explicitly reduced.
The initial state decomposition can be further simplified using the fact that it is a product state of two single photons localizing at different parts of the waveguide. Therefore the problem is reduced to decompose a single-photon state for that part of the waveguide. Imagine a single photon is localized at n adjacent time bins with the state P n i¼1 f iĉ The decomposition algorithm is as follows: (1) For the first step, calculate the density matrix for the first time-bin and diagonalize it to get its eigenvalues l ½12 a1 and eigenvectors C ½1 a1 E where a 1 is the label of different eigenvalues. The Schmidt vector l [1] is the square root of the eigenvalues, which can be truncated by keeping all those values of a 1 such that l ½1 a1 is larger than a certain threshold. The elements of tensor G [1] can be obtained by expanding the eigenvectors in the local basis {|0i, |1i}, that is, Here i 1 can take the value of 0 and 1 and a 1 only takes the values after truncation. (2) For the kth (1okrn) step, take the composite system of 1k time-bins as a whole and calculate its density matrix. Similarly, calculate the square root of the eigenvalues and then truncate to get the Schmidt vector l [k] . Take the inner product C ½k À 1 , where |i k i is the local basis {|0i, |1i} of the kth time-bin and a k À 1 , a k are the index of the truncated Schmidt vectors l [k À 1] and l [k] , respectively.
MPS update. To evolve the many-body state, a sequence of unitary operatorsÛ n are applied to update the MPS. Herê As already mentioned in the main text, the MPS used for simulation of our protocol keeps track of a waveguide of finite length equivalent to a delay time of T m /2 and also the optomechanical system. The update of the MPS chain of N time-bins is illustrated in Fig. 4. First for every n from 1 to N, we applyÛ n to the optomechanical system and the nth time-bin. Then a swap gate is used to permute the order of the two bins, so thatÛ n þ 1 can be applied in the same local way as U n . By repetitively applying these steps, the optomechanical system is moved to the last site of the MPS having interacted with all of the waveguide. Then we move the optomechanical system back by performing swap gates on the optomechanical system and the nth time-bin for each n from N to 1. Note that in this step the many-body state is not evolved in time since it stays the same under swap gates. Repeating this process again implements one repetition of the protocol. The full evolution may be over multiple repetitions, N rep .
Extraction of the entanglement entropy. The canonical MPS form provides a convenient way to extract many properties of the system, like average values during evolution, correlation functions of the quantum field and so on. In addition to the phase and fidelity that can be directly calculated from the wavefunction, we are particularly interested in the entanglement entropy of various bipartitions of the whole system. These can be directly calculated from the l [k] , the Schmidt vector at time t k , in the MPS representation. We calculate which gives us the entanglement between part of the system, containing time-bins 0otot k , and another part, with time-bins t k otoT m /2.
Relation to a closed system. To simulate the effect of many reflections by the mirror, we use a series of swap gates to move the optomechanical system back every time it reaches the last site of the MPS. This scheme seems reasonable, and we justify it here in a more rigorous way. The system as modelled is simply a closed system without losses. Here we wish to show formally the equivalence between our simulation scheme and the dynamics of the closed system. Consider a long waveguide of length L with linear dispersion, then H B ¼o 0 P 1 À 1 nĉ w nĉn where o 0 ¼ pc/L is the frequency of the fundamental mode and c is the speed of light. The summation starts from À N since we are in a rotating frame of cavity frequency. In this model the delay time of the coherent feedback is T ¼ 2L/c.
After going into the rotating frame ofĤ B , the interaction term iŝ r X nâĉ y n e ino0t Àâ yĉ n e À ino0 t ¼ ffiffiffi k pâĉ y ðtÞ Àâ yĉ ðtÞ ; where we have definedĉðtÞ¼ 1 ffiffiffi ffi 2p p P nĉn e À ino0t . It is important to note that cðtÞ¼ĉ t þ 2p=o 0 ð Þ 1 cðt þ TÞ which is a direct result of the coherent feedback. Next let us discretize time in small steps of Dt ¼ T/N, where N is a large integer and also the number of time-bins. We define in the time domain ti À Dtĉ ðsÞds, where i is an integer and t i ¼ iDt. Sinceĉ i þ N ¼ĉ i for any i, all the independent operators areĉ i ; i¼1; Á Á Á ; N f g . As before, it is straightforward to show that ½ĉ i ;ĉ w j ¼d ij ; 81 i; j N. Therefore we have N independent harmonic oscillators as time-bins.
To describe the evolution of this state, we start again with the single-step evolutionÛ Hereĉ n is the operator defined in the time domain andÛ n ¼Û n þ N . This cyclic property of these unitary matrices evolving the state means that the total evolution operator Q jÛj over many steps can be written as for t ¼ nT þ mDt where m, n are two integers and 1 m N. When m goes from 1 to N, the optomechanical system moves through the MPS and interacts with the photons. When m becomes 1 again, a series of swap gates are required to move the optomechanical system back to the first site of the MPS, so thatÛ 1 can be applied on the MPS in a local way. This is exactly the simulation scheme we have implemented as described above.
Simulation parameters. The cutoff of the singular values in the singular value decomposition calculation is set to be 10 À 4 . Since the coupling between the cavity and the waveguide is the fastest process in the rotating frame, we choose Dt ¼ 1/(2k). The temporal width of the photon is chosen to be t ¼ 1,000/k, that is, most of the photon is spread across 1,000 Â 2 ¼ 2,000 time-bins to ensure that it can fully enter the cavity. We choose the cutoff of the Fock space to be 2 for the optical cavity, 15 for the mechanical oscillator and 2 for each time bin. The total waveguide is composed of T m /(2Dt) time-bins, which is 41,888 for o m /k ¼ 1.5 Â 10 À 4 . We have verified that these paramters are sufficient for capturing the physics of the problem.
Evolution of fidelity over many bounces. The quantum optomechanical nonlinearity can cause the state of the quantum field in the waveguide to move outside of the subspace defined by the four state vectors |jki, where there are j, k ¼ 0, 1 photons in the temporal modes j and k of the waveguide. There are three effects that cause this. One is the change in the temporal mode after each bounce from the cavity. The other is leakage due to the nonlinearity. Finally there is the effect of residual entanglement with the mechanical system after a repetition of the protocol. ......

Swap
Step 1 Step 2  Step 1: for each time-bin n from 1 to N (the length of the chain), apply gateÛ n to the optomechanical system and the nth time-bin, followed by a swap gate to permute their order.
Step 2: use a series of swap gates to bring the optomechanical system back to the beginning of the matrix product state chain, so that steps 1 and 2 can be repeated again to complete one repetition of the protocol.
We are interested in the latter two since the first can be in principle reversed with linear optical components. We calculate the change in the overlaps |hjk|Ci| 2 starting from the initial state |Ci i as defined in the main text, and with |jki evolved separately with g 0 ¼ 0. The results are plotted in Fig. 5. Nominally, F ¼ 0.25 in the perfect case for every j and k, and so the infidelity is 1 À 4F.
Semiclassical calculation of fidelity. The fidelity of the CPHASE gate depends on the input state. In a semiclassical model in the bad cavity regime and with photons with sufficient temporal extent t, the fidelity of the gate for |00i, |10i and |01i input states is close to one. However, for |11i input state, the cavity frequency shift due to the interaction of the first photon reduces the probability of the second photon entering the cavity and causes a residual phonon number at the end of the protocol.
We can estimate this effect in a semiclassical model. A photon bouncing off a mechanical system that is in a coherent state |bi (instead of in its ground state) can be described by a displacement operator in mechanical phase space: The reason is that the cavity shift induced by the mechanical displacement (b þ b*) reduces the momentum kick induced by the interaction with the photon. After the interaction with first photon, the mechanical changes from being in its ground state |0i to U b ð0Þ 0 j i¼ À ir j iwhich is also a coherent state, where r ¼ 4g 0 /k. After T m /4 free evolution, the mechanical state becomes À r j i-the momentum kick induced by the photon becomes a displacement. Then the interaction with the second photon changes the mechanical state to U b ð À rÞ À r j i¼e if 1 j À rð1 þ i 1 þ r 4 Þi, where f 1 ¼ r 2 1 þ r 4 . Continuing this calculation until the end of one repetition of the protocol leads to a mechanical final state b r j i with b r ¼ r 5 1 þ r 4 ð Þ 2 þ r 4 þ ir When ro1, b r Er 5 (1 þ i). To compare the contribution of this effect to the fidelity of our gate, we numerically keep track of the evolution of the mechanical state for 11 j i input state and calculate the fidelity as well as the phase shift for each repetition of the protocol. As shown in Fig. 6, the semiclassical effect is a major part of the contribution to the full quantum results. Here the fidelity is defined as 0jb r h i j j 2 ¼e À b r j j 2 .
Fidelity change caused by another mechanical mode. The Hamiltonian of system with other mechanical modes will have other terms x(t)a w a where x(t) ¼ g p x p (t) is the frequency jitter due to a parasitic mechanical mode. We assume there is only one other mechanical modes, that during a photon interaction, x(t) is approximately constant, and model this parasitic mode classically. Then the extra phase shift can be estimated as exp À ix R t a y ðtÞaðtÞdt ¼ exp À i4x=k ð Þ . With an entangled photonic state input þ j i¼ 1 ffiffi 2 p 10 j iþ 01 j i ð Þ , after a perfect phase gate, since 10 j i and 01 j i should have the same phase, the state will remain þ j i. However, since the parasitic mechanical mode incurs a random phase shift, the final state becomes 10 j iþe iy 01 j i with y¼ 4 k x 1 þ x 3 À x 2 À x 4 ð Þ . Here x 1 (x 2 ) and It is worth noting that the agreement between the simple classical calculation and the quantum model becomes worse at larger coupling g 0 /k and that quantum model seems to lead to larger phase shifts. These are results of simulations done for g 0 /k ¼ 0.1, o m /k ¼ 1.5 Â 10 À 4 , t ¼ 1000/k and N rep ¼ 1,y,10. Note that for |00i subspace the results are always 0 and therefore not shown here.
x 3 (x 4 ) are related to the random phase shift obtained by the 10 j i 01 j i ð Þ during its two bounces on the mechanical oscillator. We model x i as independent random variables obeying normal distribution N 0; s 2 ð Þ. Therefore the distribution for y is y $ N 0; 4 4s k À Á 2 . For a given y the fidelity is FðyÞ¼ 1 þ cos y 2 , and on average the fidelity is F¼ Z dyf ðyÞFðyÞ¼ Z dyf ðyÞ 1 þ cos y 2 ð21Þ where f(y) is the probability distribution of y. Since the phase fluctuation can be estimated by s¼ gp op ffiffiffiffiffiffi kB T mp q , we arrive at an expression bounding g p and the temperature T for a given required fidelity.
Data availability. The data and code relating to this work are available on request from the corresponding author.