Advanced-Retarded Differential Equations in Quantum Photonic Systems

We propose the realization of photonic circuits whose dynamics is governed by advanced-retarded differential equations. Beyond their mathematical interest, these photonic configurations enable the implementation of quantum feedback and feedforward without requiring any intermediate measurement. We show how this protocol can be applied to implement interesting delay effects in the quantum regime, as well as in the classical limit. Our results elucidate the potential of the protocol as a promising route towards integrated quantum control systems on a chip.

In advanced-retarded (A-R) differential equations, or mixed functional differential equations, the derivative of the associated function explicitly depends on itself evaluated at different advanced-retarded values of the variable [1][2][3][4][5] . In order to solve such A-R equations either analytically or numerically, we require the knowledge of the solution history out of the domain of the equation. In many scientific disciplines A-R differential equations are used to describe phenomena containing feedback and feedforward interactions in their evolution [6][7][8] . In physics, for instance, A-R equations can be used to model dynamical systems exhibiting certain symmetry in the evolution. As a prominent examples one may mention the application of A-R equations 9 in Wheeler-Feynman absorber theory 10,11 and in the propagation of waves in discrete spatial systems 12 .
In the context of Quantum Mechanics, the implementation of feedback is more intricate than in the classical case due to the sensitivity of quantum systems to measurements. In this regard, a set of techniques has been developed for the realization of feedback-dependent systems, each of them employing different resources such as dynamical delays [13][14][15] , machine learning optimization 16 , weak measurements 17,18 , including quantum memristors 19,20 , and projective measurements for digital feedback 21 , among others. Certainly, the inclusion of feedback or memory effects in quantum dynamical systems has extended the scope of quantum protocols, and it has allowed for the study and reproduction of more complex phenomena. Therefore, devising schemes for engineering Hamiltonians that display advanced-retarded dynamics is of great relevance. Along these lines, the field of non-Markovian quantum dynamics focuses on the study of effective equations that govern the evolution of systems interacting with environments depending on previous times [22][23][24][25] . As a result, estimation of non-Markovianity sheds light on the memory content of the systems under study.
In this article we show that photonic lattices can be used to effectively tailor the dynamics of classical and quantum light fields in an advanced-retarded fashion. Our strategy is to exploit the duality between light propagation in space and time evolution 26 . Our A-R photonic approach exploits the isomorphism existing between the steady state of judiciously-designed photonic waveguide circuits and solutions of A-R differential equations. We foresee that the inherent versatility of the proposed system will make the implementation of feedback and feedforward noticeably simple, in both quantum and classical frameworks, and thus may pave the way to interesting applications in integrated quantum technologies.

Results
Advanced-Retarded Analogy. In order to introduce our protocol, we start by considering the following first order linear and non-autonomous A-R differential equation with κ + (t) = κ − (t + τ), and boundary conditions κ − (0 < t < τ) = κ + ((N − 1)τ < t < Nτ) = 0. Associating the functions x(t), x(t ± τ) with the mode amplitude of monochromatic waves traversing the j-th and j ± 1-th waveguides, a j (z) and a j±1 (z), of an array of N evanescently coupled waveguides, each supporting a single mode and having a "time" dependent propagation constant β(z), we obtain a system governed by a set of N differential equations In the quantum optics regime, the dynamics of single photons traversing this type of devices is governed by a set of Heisenberg equations that are isomorphic to Eq. (2). The only difference is that in the quantum case, the mode amplitude a j is replaced by the creation operators † a j 27 . In order to make Eq. (2) isomorphic to Eq. (1), we must impose a continuity condition a j (0) = a j−1 (τ) within the interval j ∈ [2, N − 1]. Physically, this condition implies that the mode field a j−1 , at the propagation distance of (j − 1)τ (τ is the length of the waveguides), has to be fed back into the input of the j-th waveguide. Furthermore, the finiteness of the waveguide array imposes the boundary conditions κ 1,0 = κ N,N+1 = 0. Additionally, we establish the aforementioned mapping of the independent variable t into the spatial coordinate z. Therefore, our protocol requires the implementation of waveguide lattices endowed with input-output connections as illustrated in Fig. 1. We point out that the phase introduced via the fiber propagation can be traced and disregarded given that it can be made equal for all fiber connections, being therefore a global phase.
The analogy we present here can be exploited with classical electromagnetic fields or single photons. The motivation for doing a quantum simulation is the possibility to embed the advanced-retarded dynamics in more general quantum protocols allowing the implementation of quantum feedback and feedforward. In the classical case, the initial condition a 1 (0) is implemented by continuously injecting light into the system. This condition is crucial to establish the isomorphism between the light dynamics in the waveguide array and the solution of Eq. (1). As a result, the solution of Eq. (1) is obtained in the stationary regime of our photonic system. Once the intensity is measured, the modulus square of the solution is obtained by merging the intensity evolution of each waveguide in a single variable. See Fig. 2 for a demonstration of the potentiality of our protocol, in which we analyze the setup depicted in Fig. 1. The light dynamics occurring in such an array is governed by Eq. (2).
An interesting point to highlight here is the existence of a z reversal symmetry in the simulation with respect to the central point of the evolution, z c = Nτ/2, for constant lattice parameters. This relation holds for the modulus square of the solution, aa*(z c + z) = aa*(z c − z). Consequently, after the system reaches the steady state, it simultaneously fulfills the periodic boundary condition, a(0) = a(Nτ), where we have neglected a global phase factor. This property combined with the space-time analogy opens a possible framework for the study and implementation of Closed and Open Timelike Curve gates [28][29][30][31][32][33][34] .
Non-Autonomous Equations. A more general scenario arises when considering space-dependent parameters in Eq. (2), which in the context of photonic lattices, means that the system becomes dynamic. The dynamic character is achieved with modulations of the refractive index of individual waveguides. We consider an implementable system conformed by a periodic variation β(z) = β 0 + ε cos(ωz), where β 0 is a constant, ε is the modulation amplitude, and ω stands for the modulation frequency along z. For illustration purposes, we provide a numerical calculation of this particular photonics simulator in the Supplementary Fig. 1. Notice that in the context of advanced-retarded equations the time dependence allows us to encode non-autonomous equations, which are hard to compute in general.

Systems of Equations.
We now turn our attention to demonstrate how our protocol can be extended to provide solutions of systems of A-R equations. To this end, we encode every unknown function in a waveguide array, and place all the involved arrays close to each other in such a way that light fields traversing the system can tunnel from array to array. In this manner, all functions are self-coupled and coupled to others, enriching the dynamics of the systems. In order to explain the operating principle of the protocol, we focus on the simplest case of two variables. Note that in this case, we can relate each array to a component of a qubit to be implemented in the dual-rail encoding with a single photon. As a first possibility, we can put together two lattices in which the feedback takes place as depicted in Fig. 3a. In this scenario, the light can hop to neighboring sites as well as to the sites of the adjacent array. This arrangement enables the time evolution simulation of a single qubit Hamiltonian combined with two terms corresponding to advanced and retarded couplings  A second configuration arises by mixing the connectors as shown in Fig. 3b. In this case we flip the qubit in the advanced and retarded times. Even though the equation exhibits the same structure as Eq. (3), the forward and backward Hamiltonians, H + and H − are the same matrices given in Eq. (4) multiplied to the left by σ x . Here H(t), H + (t) and H − (t) depend on the propagation constants, on the coupling between waveguides belonging to different arrays and the coupling between the guides of the same array. The vertical coupling is refereed to as q while the transverse coupling constants are represented by κ and the labels (x, y) correspond to the plane where the arrays are located. Moreover, d accounts for the diagonal coupling. In addition to the conditions for each array in Eq. (1), the following consistency and boundary conditions are fulfilled, In Supplementary Fig. 2 a numerical calculation of Eq. (3) is provided. Even though our formalism is based mainly on classical optics, one can emulate a quantum system (qubit) via the encoding of multiple A-R equations, as shown in this section. Moreover, one may also include quantum effects via the loading of the chip with genuinely quantum photonic states, e.g., two-photon or few-photon Fock states, which will introduce quantum effects as photon bunching, in combination with advanced-retarded physics. A thorough analysis of this last possibility is, however, outside of the scope of this article.
Multiple Delays. We next consider a variant of Eq. (1) with multiple delays. This configuration arises when we allow each waveguide to couple to multiple neighbors and reordering the feedback connectors. The first non-trivial example is a two-time A-R equation, given by Experimentally, the arrangement can be engineered by fabricating the waveguides in a zig-zag configuration. The resulting equation shares the structure of an A-R differential equation with additional feedback and feedforward terms. For this particular system the coupling coefficients are related as with the appropriate boundary conditions. See Supplementary Fig. 3 for a numerical simulation of Eq. (6).
Higher Order Equations. One last generalization of the A-R simulator consists in introducing complex dynamical parameters. This can be achieved by combining our feedback technique with Bloch oscillator lattices [35][36][37] . These types of arrays can be implemented by including a transverse ramping on the potential of the waveguides or by curving the waveguide arrays. Provided the evolution equations for the Bloch oscillator array, 1 . The inclusion of arbitrary complex parameters could be used to enhance the versatility of the protocol. Furthermore, a complete control of the coupling constants would allow to simulate higher order equations via systems of first order equations.
Even though the toolbox introduced here is valid for simulating diverse physical configurations, it can be generalized by an appropriate mathematical treatment. Let M(t) be the matrix containing the information about the propagation constant and couplings among the waveguides defined by =  ia t M t a t ( ) ( ) ( ), F the matrix encoding the input-output connections, and α the initial state independent of the feedback and feedforward mechanism, such that a(0) = Fa(τ) + α. Consider now that the dynamical equation is solvable in terms of the evolution operator U(t), a(t) = U(t)a(0). Our goal is to determine the complete initial condition a(0) in terms of the independent initial condition α, evolution operator U(t) and input-output matrix F. The consistency relations at a(τ) impose a set of equations that a(0) has to fulfill, a(τ) = U(τ)a(0), Notice that this relation holds for any α, allowing the input of quantum states superposed in more than one waveguide, and is also valid for different configurations of couplings U and connections F, limited by the existence of the inverse of ( − FU(τ)). Moreover, we can think of different experimental conditions, in which the connections happen at distinct evolution times τ i , Implementation Analysis. We analyze now the main limitations of the protocol regarding its interpretation and experimental realization.
Regarding the quantum equation, although the probability is not normalized during the dynamics, it is normalized in the input and output points. Therefore, initialization and measurements retrieve the correct interpretation, even if the particle undergoes forward and backward jumps on time, an effect that could be useful for the simulation of absorbing potentials 38 .
Stationary State Proximity. A natural source of errors is given by the waveguide losses damaging the quantum state in the time elapsed until the photon escapes from the chip. This time is related with the population in the solution, which is unknown before the experiment is realized. Identifying the relation between the resonances and the dynamical parameters could be helpful for estimating the population in the stationary state, and therefore the total experimental time for achieving the stationary state in the classical case. See Fig. 4 for a scheme of the error depending on the distance from the stationary state.
Experimental Errors. We have to take into account that the losses introduced by the fiber connections will break the continuity condition allowing us to simulate Eq. (1) in terms of Eq. (2). The length and propagation constant of this fibers have to be tuned so that no phase is introduced in the evolution. Additionally, the space dependence of propagation and coupling constants is limited by the experimentally realizable functions. The degrees of freedom to be considered are the writing precision for modifying the propagation constants and the spatial path of each waveguide for modifying the coupling constants.

Discussion
In conclusion, we have developed a flexible but realistic toolbox for implementing advanced-retarded differential equations in integrated quantum photonics circuits. We have shown that our analogy enables the simulation of time dependent systems of advanced-retarded equations, which in the context of quantum information can be employed to realize feedback and feedforward driven dynamics. Therefore, we consider that our work enhances the versatility of quantum simulators in the abstract mathematical direction and in terms of applications for retrospective and prospective quantum memory. We depict the decimal logarithm of the error as a function of time for three runs of the simulation with different distance with respect to the stationary state. The fact that the effective interaction between photons is null makes possible the analogy between the stationary state solution and the accumulation of solutions for an initial excitation combined until the initial population has escaped from the output port. Therefore, the distance is calculated as the norm of the population that remains in the chip. The dynamical constants of the system are equivalents to the ones in Fig. 2.