Time evolution with symmetric stochastic action

Quantum dynamical time-evolution of bosonic fields is shown to be equivalent to a stochastic trajectory in space-time, corresponding to samples of a statistical mechanical steady-state in a higher dimensional quasi-time. This is proved using the Q-function of quantum theory with time-symmetric diffusion, that is equivalent to a forward-backward stochastic process in both the directions of time. The resulting probability distribution has a positive, time-symmetric action principle and path integral, whose solution corresponds to a classical field equilibrating in an additional dimension. Comparisons are made to stochastic quantization and other higher dimensional physics proposals. Five-dimensional space-time was originally introduced by Kaluza and Klein, and is now widely proposed in cosmology and particle physics. Time-symmetric action principles for quantum fields are also related to electrodynamical absorber theory, which is known to be capable of violating a Bell inequality. We give numerical methods and examples of solutions to the resulting stochastic partial differential equations in a higher time-dimension, giving agreement with exact solutions for soluble boson field quantum dynamics. This approach may lead to useful computational techniques for quantum field theory, as the action principle is real, as well as to ontological models of physical reality.


I. INTRODUCTION
The role that time plays in quantum mechanics is a deep puzzle in physics. Quantum measurement appears to preferentially choose a particular time direction via the projection postulate. This, combined with the Copenhagen interpretation that only macroscopic measurements are real, has led to many quantum paradoxes. Here, we derive a time-symmetric, stochastic quantum action principle to help resolve these issues, extending Dirac's idea [1] of future-time boundary conditions to the quantum domain. In this approach, quantum field dynamics is shown to be equivalent to a time-symmetric stochastic equilibration in the quasi-time of a higher dimensional space, with a genuine probability. There are useful computational consequences: an action principle with a real exponent has no phase problem.
The theory uses the Q-function of quantum mechanics [2][3][4], which is the expectation value of a coherent state projector. It is a well-defined and positive distribution for any bosonic quantum density matrix, and can be generalized to include fermions [5]. The corresponding dynamical equation is of Fokker-Planck form, with a zero trace, non-positive-definite diffusion. This leads to an action principle for diffusion in positive and negative time-directions simultaneously, equivalent to a forward-backwards stochastic process. The result is timereversible and non-dissipative, explaining how quantum evolution can be inherently random yet time-symmetric.
Using stochastic bridge theory [6][7][8], the Q-function time-evolution is shown to correspond to the steady-state of a diffusion equation in an extra dimension. Thus, stochastic equilibration of a classical field in five dimen- * peterddrummond@gmail.com sions gives quantum dynamics in four dimensional spacetime. This shows that classical fields in higher dimensions can behave quantum-mechanically, including all the relevant real-time dynamics. No imaginary-time propagation is required, and the statistical description is completely probabilistic. A treatment of measurement theory is given elsewhere, showing that with this approach, a projection postulate is not essential, as only gain is needed to understand measurement [9,10].
For the fields used here to be equivalent to quantum fields, they must propagate stochastically in a negative as well as a positive time-direction. Time symmetric evolution was proposed by Tetrode in classical electrodynamics [11]. Dirac used the approach to obtain an elegant theory of classical radiation reaction [1], which was extended by Feynman and Wheeler [12]. Time-reversible methods are also studied in quantum physics [13][14][15][16], the philosophy of science [17], and used to explain Bell violations [18]. Here, we use this general approach to analyze interacting fields, thus giving time-symmetric quantum physics a strong theoretical foundation.
By comparison, the Fenyes-Nelson approach to stochastic quantum evolution [19,20], does not have a constructive interpretation [21]. The approach of stochastic quantization [22] uses imaginary time. Such methods have the drawback that analytic continuation to real time dynamics can be intractable [23,24]. The mathematical technique used here combines the Wiener-Stratonovic stochastic path integral [25,26], with Schrödinger's [6] idea of a stochastic bridge in statistical mechanics, as generalized by later workers. The resulting classical equilibration is exactly equivalent to quantum dynamics.
All quantum effects are retained in this approach, including Bell violations [10]. This is not unexpected, because quantum absorber theory, with similar timereversed propagation, also has Bell violations [15]. The focus of this paper is to understand quantum dynamics and measurement using stochastic methods. This is important both for fundamental applications to quantum measurement theory [10]. In addition, stochastic methods scale well in computation involving large systems. This may therefore help to compute exponentially complex many-body dynamics.
The Kaluza-Klein theory of electromagnetism [27][28][29], string theory [30,31] as well as the Randall-Sundrum [32] and Gogberashvili [33] approach to the hierarchy problem all use extra space-time dimensions. In the present theory the extra dimension is time-like and non-compact. Although it is not necessary to take this literally, one could ask at which coordinate in the fifth dimension is our universe? This is answerable in anthropometric terms. Just as in 'flatland' [34], the location of observers defines the extra coordinate. It is not impossible to generalize this approach to Riemannian metrics.
The Q-function is probabilistic and defined in real time. Yet it does not have a traditional stochastic interpretation, since unitary evolution can generate diffusion terms that are not positive-definite [35]. An earlier method of treating this was to double the phase-space dimension to give a positive diffusion [36]. This is usually applied to normal ordering [37], but the corresponding distribution is non-unique, and is most useful for damped systems [38] or short times [39][40][41]. With undamped systems, doubling phase-space gives sampling errors that increase with time [42,43]. Rather than using this earlier approach, here a positive diffusion is obtained through equilibration in an extra space-time dimension.
Quantum dynamical problems arise in many fields, from many-body theory to cosmology. The utility of the path integral derived here is that it is real, not imaginary [44]. Other methods exist for quantum dynamics. These include mean field theory, perturbation theory, variational approaches [45], standard phase-space methods [3] and the density matrix renormalization group [46]. Each has its own drawbacks, however. The time-symmetric techniques given here use a different approach, as well as providing a model for a quantum ontology.
To demonstrate these results, we introduce a general number conserving quartic bosonic quantum field Hamiltonian. The corresponding Q-function dynamics satisfies a Fokker-Planck equation with zero trace diffusion. This leads directly to a time-symmetric action principle. The corresponding probabilistic path integral has a solution obtained through diffusion in a higher dimension. Elementary examples and numerical solutions are obtained. We compare results with exactly soluble cases.
The content of this paper is as follows. Section II summarizes properties of Q-functions, and proves that they have a traceless diffusion for number conserving bosonic quantum field theories. Section III derives the action principle. Section IV treats extra dimensions, and shows how the classical limit is regained. Section V gives examples and numerical results. Finally, section VI summarizes the paper.

II. Q-FUNCTIONS
Phase-space representations in quantum mechanics allow efficient treatment of large systems via probabilistic sampling [47]. These methods are very general. They are related to coherent states [37] and Lie group theory [48], which introduces a continuous set of parameters in quantum mechanics. Results for bosonic fields are summarized in this section. The Q-function method [2] can also be used for spins [49,50] and fermions [5] as well as for bosons, with modifications. These cases are not treated in detail here, for length reasons.

A. General definition of a Q-function
A general abstract definition of a Q-function [5] is: whereρ (t) is the quantum density matrix andΛ (λ) is a positive-definite operator basis, and λ is a point in the phase-space. This must give an expansion of the Hilbert space identity operatorÎ, so that, given an integration measure dλ,Î = Λ (λ) dλ . The basis is not orthogonal, and it is generally essential to employ non-orthogonal bases and Lie groups in order to obtain differential and integral identities. Provided T r [ρ (t)] = 1, the Q-function is positive and normalized to unity: It therefore satisfies the requirements of probability. Quantum expectations Ô Q of ordered observablesÔ are identical to classical probabilistic averages O C -including corrections for operator re-ordering if necessary -so that: Here, Q indicates a quantum expectation value, C is a classical phase-space probabilistic average, and timearguments are implicit. The basis functionΛ does not project the eigenstates of a Hermitian operator, and therefore the quantum dynamical equations differ from those for orthogonal eigenstates.
The examples treated here use the Q-function for a complex N −component bosonic fieldψ (r). This is defined with an n d -dimensional space-time coordinate r, where r = r 1 , . . . r n d = (r, t). Quantum fields are expanded using M annihilation and creation operatorŝ a i ,â † i for M/N spatial modes. These describe excitations localized on a spatial lattice, or single-particle eigenmodes. The indices i include the N internal degrees of freedom like spin quantum numbers and/or different particle species.
On Fourier transforming to position space, Q [ψ] in field space is a functional of the complex field amplitudes ψ (x). Results can be calculated in either field or mode notation. Either approach is equivalent in terms of the resulting dynamics. From now on we focus on the mode expansion method, noting that there are equivalent formulations using functional integrals.

B. Observables
The transition probability or expectation of any observableσ is obtained by expandingρ in a generalized P-representation, P (α, β). This always exists [52], so that for any quantum density matrixρ, (2.5) HereΛ p (α, β) is an off-diagonal coherent projector, and dα, dβ are each M dimensional complex integration measures, so that if α = α x + iα y , then dα = d M α x d M α y . The existence proof [52] shows that there is a canonical probability distribution P (α, β) given by: (2.7) We now show that this leads to a general operator correspondence function forσ in the form of: To prove this, we use the expansion ofρ in Eq (2.5), which gives that: Expanding this using the canonical expansion, Eq (2.7), the c-number function corresponding toσ is therefore O σ (α), where: (2.10) As a simple example, particle numbers in the bosonic case are given by introducing the equivalent c-number function n (α) ≡ |α| 2 − 1, so that the quantum and classical averages agree: This is a special case of the more general identity given above. As another example, when expanded in mode operators, an n − th order anti-normally ordered moment is: The operator moments can be of any order.
Similar techniques are available for fermions [5] and spins [49,50], so this approach is not restricted to bosonic fields. As first emphasized in Dirac's early review paper [53], one can calculate any observable average from a classical looking distribution, provided the observable is re-expressed in terms of a suitable operator ordering. In this case, it is the anti-normal ordering of ladder operators that is utilized, and the resulting distribution is always positive.

C. Exact results and identities
Exact analytic solutions for the Q-function are known for a number of special cases, including all gaussian states. For a noninteracting multi-mode vacuum state, |Ψ = |0 , and more generally for a coherent state |Ψ = |α 0 c , where α 0 ≡ 0 in the vacuum state, one obtains This has a well known interpretation [54,55]. If one makes a simultaneous measurement of two orthogonal quadratures, which is possible using a beam-splitter, then Q (α) is the probability of a simultaneous measurement of quadratures α x and α y , where α = α x + iα y . This is also the result of an amplified measurement [56]. Similarly, any number state |Ψ = |n has a simple representation as: (2.14) A free-particle thermal state with mean particle number n th is given by: There are several mathematical properties that make this expansion a very interesting approach to quantum dynamics. We first introduce a shorthand notation for differential operators, ∂ n ≡ ∂/∂α n . There are the following operator correspondences [4,37,57,58] written in terms of mode creation and annihilation operators: Q-function evolution equations are obtained by using the operator identities to change Hilbert space operators acting onρ to differential operators actingΛ, and hence on Q α .
There is a product rule for operator identities. The full set of 2M mode operators are written with a superscript asâ µ for µ = 1, . . . 2M , whereâ j =â j ,â j+M =â † j , and L µ denotes the corresponding differential term. The general identities can be written: (2.17) To obtain operator product identities, one uses the fact that the mode operators commute with the c-number terms, so that the operator closest to the kernelΛ always generates a differential term that is furthest from Λ:â

D. Quantum field dynamics
To understand time-evolution, we consider an arbitrary time-dependent multi-mode Hamiltonian with quartic, cubic and quadratic terms and a generalized number conservation law, typical of many common quantum field theories. This generic quartic Hamiltonian is re-expressed by expanding all fields with mode operators. Here we choose to use an antinormally ordered form, for simplicity in applying the operator identities, which therefore gives us that: (2.19) For notational convenience, to combine all the terms in one sum, we include summations over i = 0, . . . M , and defineâ 0 = 1. While formally quartic, this includes linear, quadratic and cubic terms as well, through the terms that includeâ 0 . This is the most general quartic number preserving Hamiltonian, which has no more than quadratic terms in either creation and annihilation operators. The time argument of g ijkl (t) is understood if it is not always written explicitly. Without loss of generality, we assume a permutation symmetry with (2.20) AsĤ must be hermitian, We assume a momentum cutoff, so that any renormalization required is carried out through the use of cutoff dependent coupling constants. Cubic terms of form a 0âjâ † kâ † l =â jâ † kâ † l are included, and these describe parametric couplings that have a generalized form of number conservation. From the Schrödinger equation, , so the dynamical evolution of the Q-function for unitary evolution is immediately given by: After implementing the mappings given above, we obtain: and defining L HΛ (α) =ĤΛ (α) as the mappings of operators on the left, andL HΛ (α) =Λ (α)Ĥ for operators on the right, we sequentially use the identities of Eq (2.16) and (2.18) to expand this as: Here we define α 0 = 1, ∂ 0 = 0, so that all cases are included. This gives the general differential equation, Similar results hold for the functional approach, but the mode expansion approach is used here for its greater simplicity.
The main thrust of the present paper is to treat unitary systems, as described by the equations above. One can include decoherence and reservoirs by adding them to the Hamiltonian. Although these can also be treated with a master equation, we assume here that any reservoirs are simply included in the dynamical equations.
Next we define an extended vector α µ , and corresponding derivatives ∂ µ for µ = 1, . . . 2M , where α j = α j , α j+M = α † j , which includes amplitudes and conjugates. Using an implicit Einstein summation convention over µ, ν = 1, . . . 2M , we note that constant terms cancel, and: From using Eq (2.21) and (2.16), the diffusion term for 1 ≤ k, l ≤ M is: Letting µ , ν ≡ µ−M, ν −M , we find that D µν α = D µ ν * α , and for unitary evolution there are no cross-terms D µν α . Similarly, using permutation symmetry, the drift term is: 29) and the conjugate drift for µ > M is A µ = A µ * α . Generally, the second-order coefficient D µν α (α) depends on the phase-space location, except for a purely quadratic Hamiltonian. In this case, the diffusion is either zero or constant in phase-space.

E. Traceless diffusion and time-reversibility
We now prove that for unitary quantum evolution, the diffusion matrix is divided into two parts, one positive definite and one negative definite, corresponding to diffusion in the forward and backward time directions respectively. We first show that the corresponding Q-function time-evolution follows a Fokker-Planck equation with a traceless diffusion matrix. That is, the equation has an equal weight of positive and negative diagonal diffusion terms.
This was previously demonstrated in studies of a different type of Hamiltonian, for thermalization in Q-function dynamical equations for spin systems [59]. The result is also true for Bose and Fermi quantum fields, and is generic to second-order unitary Q-function evolution equations. In this paper we show it for Bose fields only. The proof in the Fermi case will be given elsewhere. To map Hilbert space time-evolution to phase-space time evolution, and to prove that the resulting diffusion term is traceless, we apply the operator identities.
Using the results of Eq (2.28), terms with non-zero i and k indices generate second-order derivative terms which combine to give the diffusion matrix. The diagonal second order terms of interest are obtained when two derivatives both act on the same mode. If all the Hamiltonian terms have k = l = 0 the diffusion is constant, but otherwise it depends on the phase-space position α.
The k −th diagonal diffusion term in complex variables ∂ k comes from from identities involvingĤ ijklΛ with 0 < k = l ≤ M , which therefore is: The phase term η(α) depends on the coupling and amplitudes. This diagonal term is accompanied by the hermitian conjugate term derived from the reverse ordering, of formΛĤ ijkl , so that L α is real overall. These conjugate terms have derivatives ∂ * j , which give, on defining This allows the introduction of real quadrature variables X j , defined such that for µ = j ≤ M : Hence, the derivative terms become: If we define X j+M = Y j , we have an extended 2M dimensional real vector, which is written with a superscript as X µ . After making this transformation, the diagonal diffusion term in real variables is: (2.34) Here k = k + M , and as a result, on summing the diagonal terms, the diffusion matrix with these variables is traceless, i.e., T r [D] = 0. This is different to classical Fokker-Planck theory, where diffusion matrices are positive-definite [60,61]. The Q-function diffusion matrix is not positive-definite, yet the distribution remains positive, from its construction as a positive-definite observable. Given this analysis, the traceless property applies to a general class of quadratic, cubic and quartic Hamiltonians. There can also be variables with zero diffusion, which are deterministic and are also traceless.
Traceless diffusion is preserved under both rescaling and orthogonal rotations: φ = OX, of the real quadrature coordinates. Since the diffusion matrix of a real Fokker-Planck equation is real and symmetric, it can always be transformed into a diagonal form D µ (φ) in the new variables φ, using orthogonal rotations. As a result, the transformed phase-space variables can be classified into two groups, having either positive or negative diffusion, with the equation: We choose the orthogonal rotation so that it results in a traceless diagonal diffusion with D µ ≥ 0 for µ ≤ M and D µ ≤ 0 for for µ > M . The 2M -dimensional phasespace coordinate is written with a superscript as φ µ for µ = 1, . . . 2M , and we use repeated greek indices to indicate implicit summation over 1, . . . 2M . This generates a characteristic structure which is universal for unitary evolution with all Hamiltonians of this form.
The phase-space vector φ is subdivided into two complementary pairs so that φ = (x, y), where the x variables have a positive semi-definite diffusion, and the y variables have a negative semi-definite diffusion. This universal characteristic structure is clearly not the usual positive-definite diffusion found in classical diffusion.
One can obtain a positive-definite diffusion in a phasespace, via an alternative approach of doubling the phasespace dimension, which is especially useful for open systems in the positive P-function representation [36,62]. This is also applicable to unitary evolution [39,40], but gives sampling errors that increase with time [63][64][65]. Such methods have successfully treated soliton quantum squeezing in photonics [41,66], and quantum soliton dissociation [67] in BEC dynamics [68]. Yet on long timescales sampling errors increase, because the distribution becomes less compact.
The method used here is to obtain an algorithm for the Q-function, rather than the P-function. Since the Qfunction is unique, sampling-error growth in time is minimized. However, a different approach to simulation is necessary, via positive-definite diffusion in a higher spacetime dimension, without changing the phase-space itself, as explained below.

F. Density-density coupling
If the Hamiltonian has only quadratic terms, the diffusion terms are either zero or constant in phase-space, as pointed out above. We now treat an alternative approach that leads to constant diffusion for the most common form of nonlinear coupling, namely density-density coupling. The result is a different definition of the phasespace variable φ, which has constant diffusion independent of φ, as well as being traceless and diagonal. This type of physics is found in the Bose-Hubbard model, and many other bosonic quantum field theories [69,70].
On a lattice we consider a quartic Hamiltonian of form: Using the identities of Eq (2.16) again, the second-order derivative terms in α are: In this case one may define a mapping, θ j = λ ln(α j ), where λ is a scaling factor, so that in the new variables the diffusion matrix D θ ij is constant, where: This transformation simplifies the derivation of the Fokker-Planck path integral. Path integrals for spacedependent diffusion exist [71], but are more complex. For the analysis given here, it is simpler to transform to the case of constant diffusion, although it is possible to obtain a path integral without doing this. If the diffusion is constant, as in a quadratic Hamiltonians, this step is unnecessary.
We now check that the traceless property persists after making this variable change to logarithmic variables. The Q-function is mapped to a set of constant diffusion, complex phase-space variables θ, as shown in Fig (1), which satisfy an equation of form: To prove the traceless property, we make a second mapping to a real quadrature vector, φ = [φ 1 , . . . φ 2M ], described by the linear transformation φ = T θ. In this constant diffusion space, there are diagonal second derivative terms together with conjugate terms such that D θ jj = e −2iηj D θ jj . The corresponding real variables are defined in this case as: This clearly results in traceless diffusion. For quadratic cases with no logarithmic transformation, one may simply define (x, y) as in the notation of the previous subsection. For a one-mode case, the mapping transformation matrix to real variables is The result is a transformed Q-function, Q = Q θ |δθ/δφ|, which evolves according to real differential equation. Introducing ∂ µ ≡ ∂/∂φ µ , for µ = 1, . . . 2M , the time-evolution equation with a diagonal, constant diffusion is: The transformed diffusion matrix is traceless as previously, so that: This means that the diagonal diffusion matrix can be subdivided into positive and negative constant diffusion parts. The phase-space probability Q is positive, yet since the overall diffusion term D is not positive definite, so this is not a stochastic process of the traditional, forward-time, type [72]. The form of Eq (2.42) means that probability is conserved at all times provided boundary terms vanish, i.e., The Q-function has very unusual properties. It is probabilistic, and obeys a generalized Fokker-Planck equation. Yet it describes a reversible process. Positive distribution functions in statistical mechanics commonly follow an equation which is irreversible, owing to couplings to an external reservoir. Despite this, the Q-function is a phase-space distribution that is positive, and has probabilistic sampling. This implies that it can be treated in a similar way to a classical probability distribution, with some modifications.

III. TIME-SYMMETRIC DIFFUSION
In the previous section, we showed that the Q-function for unitary evolution can be transformed to have a real differential equation with a traceless diffusion term. The consequence is that it has a time-symmetric diffusion which is not positive-definite. As a result, forward-time sampling via a stochastic differential equation in time is not possible. In the present section we obtain local propagators for the traceless Fokker-Planck equation. This will be used to derive an action principle and path integral. The result leads to stochastic evolution in spacetime with time-reversal symmetry.

A. Time-symmetric Green's functions
The universal property of traceless diffusion means that the real phase space of φ is generally divisible into two M -dimensional sub-vectors, so that φ = [x, y]. The x fields will be called positive-time fields, with indices in the set T + , while fields y are called negative-time fields, with indices in the set T − . The fields x are those with D µ ≥ 0, while fields y have D µ ≤ 0. These have the physical interpretation of complementary variables.
One usually solves for Q (φ, t) at t given an initial distribution Q 0 (φ, t 0 ). However, the lack of a positivedefinite diffusion in y means that one cannot use standard Green's functions to propagate Q forward in time, without requiring singular functions. Instead, the approach taken here is to solve for Q (φ, t) given an initial We therefore consider a general dynamical problem amenable to solution using time-symmetric methods. Suppose that at time t 0 the Q-function has a known initial marginal distribution for x, such that and that the final marginal distribution for y is known for a later time t f such that To solve for Q (φ, t), we define a time-symmetric , and whose negative time components end at y f = y(t f ). We abbreviate this as [φ, t |x 0 , y f ]. This is a function of φ that satisfies the generalized Fokker-Planck equation (2.42), with initial and final marginal conditions: The form of Eq (2.42) means that, using partial integration, probabilities are conserved both for Q (φ, t) and for the Green's functions, provided boundary terms vanish. We now prove three important results.

Normalization
Like the Q-function, time-symmetric Green's functions are normalized. This is because the Q-function dynamical equation conserves probability, and clearly from Eq (3.3) the [φ, t |x 0 , y f ] functions are normalized both initially and finally. As a result, for all times,

Solution theorem
If a symmetric Green's function exists for arbitrary x 0 , y f , then the solution for Q (φ, t) can be obtained by integration over x 0 and y f , so that: To prove this, note that Q (φ, t) as defined above must satisfy the Q-function differential equation, (2.42), since [φ, t |x 0 , y f ] does. One can verify through direct integration, together with the fact that the marginals P x and P y are normalized to unity from Eq (2.2), that this solution for Q also satisfies the required marginal probability conditions, (3.1) and (3.2).

Factorization theorem
Given a time-evolution equation of form (2.42), then in the limit of a short time interval ∆t = t f − t 0 , the timesymmetric Green's function factorizes into a product of forward time and backward Green's functions. In greater detail, this factorization property of the time-symmetric propagator is as follows: Definingφ ≡ [x 0 , y f ] as a vector of initial and final fields, and short time propagators x, t φ and y, t φ for the x and y fields, then φ, t φ factorizes over short time intervals as We now prove this and obtain the explicit form of x, t φ and y, t φ from the generalized Fokker-Planck equation. The relevant time-evolution equation is (2.42). The diagonal diffusion means that the Fokker-Planck equation has forward and backward parts, so the differential operator can be written Here, L + (L − ) only includes derivatives of x, (y) respectively, together with drift terms that are functions of φ.
Each is a positive-definite Fokker-Planck operator. On defining d µ = |D µ | these differential operators are: is differentiable and smooth, the drift and diffusion terms can be approximated by their initial values and times at φ [4,8,61], so that The local differential equation then has the form: Provided that boundary requirements are satisfied, this is solved by setting f (x, y) = f x (x) f y (y), where: To satisfy the boundary conditions we must impose the initial condition on , and the final condition on y that f y (y, t f ) = δ (y − y f ), while noting that the form of L ± will maintain normalization over the interval. This can be verified more rigorously by expanding A µ (φ) in a Taylor series to first order in ∆φ = φ−φ, then solving and taking the limit of ∆t → 0. Each of the differential operators corresponds to a weighted diffusion time-evolution. In the limit of a short time interval, provided each drift term is evaluated at its initial value in x (y) respectively, the time-symmetric Green's function factorizes into a product of forward time and backward time terms, where: The time interval ∆t + is in the forward time direction, measured from the start of the interval, while ∆t − is a time interval in the backward time direction, measured from the end of the interval. The normalization term is the standard one for a normalized solution to the diffusion equation [61,71]. Because of the symmetry of the diffusion matrix it is the same function in either direction, in the limit of a small time-interval: This behavior is indicated graphically in Figure (2).

B. Discrete trajectories
Consider a phase-space trajectory discretized for times t k = t 0 + ∆tdt, with k = 1, . . . n. We wish to con- For n = 1, we have an initial and final constraint, as in the Green's function boundary conditions of Eq (3.3).
One can obtain the respective probabilities of transition from x 0 → x 1 and y 1 → y 0 , over a short time interval, using the results of the last section.
Since the Q-function is a probability, these are the marginal probabilities P x (x 1 , t 1 ) and P y (y 0 , t 0 ). On integrating the factorized Green's function solution over the conjugate variables, one obtains: (3.14) From now on we drop the time arguments, taking [x 1 |x 0 , y 1 ] ≡ [x 1 , t 1 |x 0 , y 1 ], as the time intervals are equal. We obtain that the joint probability of transitions x 0 → x 1 and y 1 → y 0 both occurring, since they independent events on a short time interval, is: We can write P (y 0 , x 1 |x 0 , y 1 ) = P ([φ 0 , φ 1 ] |x 0 , y 1 ). These are equivalent, because the probability of observing x 1 given (y 1 , x 0 ) is the same as the probability of observing φ 1 , since the event at y 1 is included in the conditioning. This situation is depicted in Fig (2).
These probabilities can be extended to multiple events, as shown in Fig (3) for n = 4. Following the chain rule of probability for conditional events, the probabilities of transition from x 0 → x 1 → x 2 and y 2 → y 1 → y 0 , is therefore: This shows that the probability for a final x 2 and initial y 0 is conditioned on both x 1 and y 2 . The result for the whole trajectory is obtained by extending the argument given above recursively, since the probability for a final x k and initial y 0 is conditioned on both x k−1 and y k , so that Applying the argument n times, and definingφ k ≡ (x k , y k+1 ), this implies that: This probability can be written as a general timesymmetric stochastic action, namely: where the action is given by a sum over all the propagation weight terms: with the time-symmetric Lagrangian defined in the generalized Ito sense [72], such that the drift and diffusion terms are evaluated at the start of every step, depending on the time propagation direction for the variable in question: (3.21) To explain this limiting procedure more precisely, we use the following definitions to take the limits: To understand the physical meaning, we refer back to Fig (3), which shows four discrete segments of propagation. At the k-th step, for t = t k , the value of x k is constrained, while the value of y k is only known probabilistically. Similarly, at t = t k+1 , the value of y k+1 is constrained, but the value of x k+1 is only known probabilistically, with some probability calculated from the Green's function. As a result, the normalization of the probability is that: where the normalization has an integration measure over all coordinates except the two fixed vales of x 0 and y f : We wish to obtain a total Green's function, [φ, t |x 0 , y f ], from the trajectory probability. In the limit of ∆t → 0, and n → ∞ , we take t = t K , i.e, the target time for the Q-function can be taken as corresponding to one of the n discrete times. This is always possible to any desired accuracy, in the relevant limit. The time-symmetric Green's function over a finite interval is therefore constructed as: This is quite general, but for the remainder of the paper we treat simpler cases whose diffusion is independent of the phase-space coordinate.

C. Central difference action principle
The time-evolution equation as a forward diffusion problem for x is a well-defined initial-value problem in the positive time coordinate t + = t with a specified distribution at t + = t 0 . Similarly, the backwards diffusion problem for y is a well-defined initial-value problem in the negative time coordinate t − = t 0 + t 1 − t + , with a specified distribution at t − = t 0 .
Over every small time-interval, there is a positivedefinite Fokker-Planck differential operator in each timedirection. Each of these differential operators acts on the variables with a diffusion in the chosen time direction. Such diffusion equations have well-known path-integral solutions, from work by Wiener [25] and deWitt [73].
The important difference is that this action is timesymmetric, propagating in both time directions with complementary variables. However, there is a subtlety here. In the previous subsection, the coefficients are evaluated at the start of each interval, atφ k = (x k , y k+1 ). Stratonovich [26,74], Graham [71] and others have shown that there are other action formulae evaluated at the center of each interval, which are covariant under phasespace frame transformations. While these can be applied in all cases, the derivation is simplest if d µ = d µ (t) is independent of φ, which is the case that is derived here. A more complex result holds for general diffusion, involving the curvature tensor in phase-space with a metric equal to the diffusion matrix [71], which is outside the scope of this article.
As explained in the previous subsection, each of these differential operators corresponds to diffusive timeevolution one of the two time directions. For the case of constant diffusion there is a central-difference Green's function in the limit of a short time interval: The functions L x,y are the Fokker-Planck Lagrangian differential operators in the relevant direction. These are defined as in Eq (3.12), except with the drift term A µ evaluated at the midpoint of the time-step, so that : This adds a correction to the normalization, depending on the divergence of the drift A, resulting in an additional exponential weighting term [26,71]. With the identification that ∂φ µ /∂t ± ≡ lim dt→0 [φ µ (t ± ∆t) − φ µ (t)] /∆t, and taking the limit of ∆t → 0, the two central Lagrangians can be written as: (3.28) Because A µ ≡ A µ (x, y), there is a physical behavior that does not occur in standard diffusion. The drift in x can depend on the field y in the backward time direction, and vice-versa, as shown in Fig (3). This cross-coupling that leads to nontrivial quantum dynamics. This is not important for short times, since each infinitesimal propagator is defined relative to an initial overall φ , but it does modify the structure of long-time propagation.
The general, time-symmetric stochastic action with central differences, when the diffusion is not constant, is also obtainable. However, this involves Riemannian curvature terms in phase-space, and is outside the scope of the present article.

D. Time-symmetric stochastic action
This method will now be employed to transform interacting quantum field theory into another form, focusing on the constant diffusion case for simplicity. The probability density for quantum time-evolution is given by a path integral over a real Lagrangian, where in each small time interval the propagators factorize. Due to timereversal symmetry of the propagator, these equations can be solved using path integrals over both the propagators.
The path then no longer has to be over an infinitesimal distance in time, and the total propagators will not factorize. This is a type of stochastic bridge [6][7][8], which acts in two time directions simultaneously. Hence, to write the bridge in a unified form, with time integration in the positive time direction only, we define a combined, central difference Lagrangian as: (3.29) so that the action integral can be written in the positive time direction for t 0 < t < t 1 , with a total Lagrangian of: Here the total potential, V , includes contributions of opposite sign from the positive and negative fields, so that: This defines the total probability for an n-step open stochastic bidirectional bridge, with constant diffusion, central difference evaluation of the action, and fixed intermediate points (3.32) On integrating over the intermediate points, with drift terms defined at the center of each step in phase-space, this can be written in a notation analogous to a quantummechanical transition amplitude in a Feynman path integral. One obtains the transition probability: where t n = t 0 + n∆t, and here we only integrate over the intermediate phase-space points: The paths φ (t) are defined so that x (t 0 ) = x 0 and y (t n ) = y f are both constrained, at the initial and final times respectively. The definition of [x f , y 0 |x 0 , y f ] is the probability of both arriving at x f and starting from y 0 , given that the initial value of x is x 0 , and the final value of y is y f . Although formally similar, this propagator has a different meaning to the quantum transition amplitude, because, as it is a probability, it is always positive valued.

IV. EXTRA DIMENSIONS
Many techniques exist for evaluating real pathintegrals, both numerical and analytic. There is a formal analogy between the form given above and the expression for a Euclidean path integral of a polymer, or a charged particle in a magnetic field. Here we obtain an extradimensional technique for the probabilistic evaluation of the path integral, due to its simplicity and interesting physical interpretation. Although other methods exist, they are not investigated in this paper in the interests of brevity.

A. Equilibration in higher dimensions
Direct solutions using stochastic equations in real space-time are not feasible, but there are other methods. To make use of the real path-integral, one needs to probabilistically sample the entire space-time path, since each part of the path depends in general on other parts. To achieve this, we add an additional 'virtual' time dimension, τ . This is used in the related statistical problem of stochastic bridges, for computing a stochastic trajectory that is constrained by a future boundary condition [6-8, 75, 76].
This extra-dimensional distribution, G ([φ], τ ), is defined so that the probability tends asymptotically for large τ to the required solution: The solution is such that φ (t) is constrained so that x (t 0 ) = x 0 , and y (t f ) = y f . It has been shown in work on stochastic bridges [7] that sampling using a stochastic partial differential equation (SPDE) can be applied to cases where one of the boundary conditions is free. To define an SPDE the other boundary condition on x is specified so thatẋ (t f ) = A x (φ (t f )), with a boundary condition forẏ so thatẏ (t 0 ) = A y (φ (t 0 )). This is consistent with the open boundary conditions of the path integral in real time, since in the limit of ∆t → 0 , the path integral weight implies that one must haveẋ for almost all paths. The effect of the additional constraint vanishes as ∆t → 0, as it contributes a negligible change to the entire path integral. This condition is necessary in order to have a well-defined partial differential equation in higher dimensions.
Extra-dimensional equilibration is not used for conventional SDE sampling, as direct evolution is more efficient. However, we will show that SPDE sampling is applicable to time-symmetric propagation, where direct sampling is not possible. In this section, a simplification is made by rescaling the variables to make the diffusion d µ (t) independent of time and index, i.e., d µ (t) = d. We also assume that there is no explicit time-dependence in the Hamiltonian. The general solution is given in the Appendix.
The SPDE is obtained as follows. Firstly suppose that G ([φ], τ ) satisfies a functional partial differential equa-tion of: G .
(4.2) In order that the asymptotic result agrees with the desired expression (3.32) for G, it follows from functional differentiation of Eq (3.33), that we must define A[φ] so that: This a variational calculus problem, with one boundary fixed, and the other free. Variations vanish at the time boundaries where φ is fixed. At the free boundaries we specify thatφ = A, as explained above. At the integration boundaries both types of boundary terms are zero because they occur in terms that vanish provided that As a result, there are two type of natural boundary terms that allow partial integration to obtain Euler-Lagrange equations. Either one can set ∆φ µ = 0 to give a fixed Dirichlet boundary term, or else one can setφ µ = A µ , to give an open Neumann boundary term. This allows one to obtain Euler-Lagrange type equations with an extradimensional drift defined as: where v µ = φµ − A µ /d. The functional Fokker-Planck equation given above is then equivalent to a stochastic partial differential equation (SPDE): where the stochastic term ζ is a real deltacorrelated Gaussian noise such that

B. Coefficients
Introducing first and second derivatives,φ ≡ ∂φ/∂t andφ ≡ ∂ 2 φ/∂t 2 , there is an expansion for the higherdimensional drift term A in terms of the field timederivatives: Here, c is a circulation matrix that only exists when the usual potential conditions on the drift are not satisfied [8], while a is a pure drift without derivatives: The function U is an effective potential, which acts to generate an effective force on the trajectories: The final stochastic partial differential equation that φ must satisfy is then: The final result is a classical field equation in an extra space-time dimension with an additional noise term. It has a steady-state that is equivalent to a full quantum evolution equation, and is identical to classical evolution in real time in the zero-noise, classical limit, as shown in the next subsection.
The equations can be treated with standard techniques for stochastic partial differential equations [77], except that the equations have n d + 1 dimensions in a manifold with n d space-time dimensions. The simplest case, for a single mode, has n d +1 = 2 dimensions. In computational implementations, one can speed up convergence to the steady-state using Monte-Carlo acceleration [78].

C. Classical limit
The classical limit is for d → 0. In this limit the higher-dimensional equations are noise-free and diffusive. Including a circulation term in case potential conditions are not satisfied, one has: Substituting the classical trajectory solution,φ ν = A ν , one sees immediately that for classical trajectories, However, on this trajectory, the second derivative term simplifies to:φ 13) and therefore one obtains: This extra-dimensional equation therefore has an exact steady state solution corresponding to the integrated classical field evolution in real time, namely: Both the initial and final boundary term equations are satisfied provided one chooses x (t 0 ) = x 0 and y (t f ) = y f , if these are compatible, that is, if the dynamical equations have a solution. If one uses these equations to solve for y(t 0 ), the solution can be rewritten in a more conventional form of a classical solution with initial conditions: (4.16) The importance of imposing future-time boundary conditions in classical field problems like radiationreaction has long been recognized in electrodynamics, including work by Dirac [1], as well as Wheeler and Feynman [12]. In such theories various field components typically require future-time restrictions on their dynamics. Hence the fact that such future-time boundaries arise in the classical limit found here should not be very surprising.
Dirac [1] described his electron acceleration result in as "the most beautiful feature of the theory". He explains: "We now have a striking departure from the usual ideas of mechanics. We must obtain solutions of our equations of motion for which the initial position and velocity of the electron are prescribed, together with its final acceleration, instead of solutions with all the initial conditions prescribed." If Dirac's type of dynamical restriction is compared with the classical limit obtained here, there are are clear similarities. There is a dynamical condition that is required to derive the correct classical time evolution, obtained from a restriction on the future boundaries of the radiation field. It is a striking, even "beautiful" feature of the present approach that this type of classical, future boundary condition arises naturally from taking the relevant zero-noise, classical limit of our equations.

D. Time-symmetric stochastic differential equation
The path integrals correspond to a functional integral over stochastic paths. Hence the trajectories can be written in an alternative, intuitive form after probabilistic sampling, as a time-symmetric stochastic differential equation [79], with: (4.17) Here the two fields are propagated in the positive and negative time directions respectively, sometimes called "forward-backward" equations, while the noise terms dw are correlated over short times so that, over a small interval dt: The compelling feature of these equations is that they unify two important features: time reversibility and randomness. These types of equations also occur in stochastic control theory, and have an extensive mathematical literature proving their existence and other properties [79]. However, while they provide an insight into the structure of the stochastic equations, they cannot be readily solved using conventional algorithms for stochastic differential equations. This can be recognized by attempting to write the equations as forward time stochastic differential equations. We defineȳ (t) as a time-reversed copy of y (t), i.e., let t − = t 0 + t 1 − t, and (4.19) The stochastic differential equation that results, treating each argument as the same time t, is: Here, x (t 0 ) = x 0 andȳ (t 0 ) = y f are now both "initial" conditions, but the y coordinate is replaced byȳ instead. In other words, we can regard the stochastic differential equations as having a forward-time stochastic propagation if the drift terms include complementary fields defined at different times. However, non-locality in time prevents one from using standard, local-time algorithms for solving even these more conventional looking equations as ordinary stochastic differential equations. This behavior is not surprising, physically. If these fields had local drift terms, they would be causal, local theories that satisfy Bell's theorem, and do not correspond to quantum theory. It is possible that analytic equations like this can be used to develop a stochastic perturbation theory [80] for quantum fields. Since forward-backward stochastic equations occur in other areas as well, such techniques may have wider applicability.

E. Numerical methods
A variety of numerical techniques can be used to implement path integrals with a time-symmetric action. In this paper we solve the equivalent higher-dimensional partial stochastic differential equation with a finite difference implementation. This permits Neumann, Dirichlet and other boundary conditions to be imposed. We also explain strategies for dealing with future time boundaries, which is the most obvious practical issue with this approach. The component x propagates in the positive time direction as a random Wiener process. The expected variance for τ → ∞ is x 2 = 1 + t, with x (t, 0) = 0, x (0, τ ) = v, and v 2 = 1. Fluctuations are sampling errors due to a finite number of 10000 trajectories. Variance error bars due to sampling errors were estimated as ±2.5%, in good agreement with the difference between exact and simulated variance. A semiimplicit finite difference method [81,82] was used to integrate the equations, with step-sizes of ∆τ = 0.0002 and ∆t = 0.03. Errors from the finite step-size in τ were negligible.

SPDE integration
First, we demonstrate convergence of the higher dimensional method, using a central difference implicit method that iterates to obtain convergence at each step, including an iteration of the boundary conditions. The method is similar to a central difference method described elsewhere [77,81]. A simple finite difference implementation of the Laplacian is used to implement nonperiodic time boundaries.
In order to demonstrate convergence, Fig (4) gives the computed numerical variance in an exactly soluble example of a stochastic differential equation with no drift term. We treat one variable and C = a = 0, using a public-domain SPDE solver [83] with a random Gaussian initial condition of x(t = 0) = v where v 2 = 1, so that: This is a case of pure diffusion, where one expects the final equilibrium solution as τ → ∞ to be x 2 = 1 + t. From Eq (4.10), the corresponding higher-dimensional stochastic process has boundary conditions of x(t = 0) = v andẋ(t = t f ) = 0, while satisfying a stochastic partial differential equation: From the numerical results in Fig (4), the expected variance is reached uniformly in real time t after pseudotime τ ∼ 2.5, to an excellent approximation, reaching x 2 = 1.95 ± 0.05 at t = t f = 1 and τ = 5.
For the examples given here, our focus is on accuracy, not numerical efficiency. The purpose of these examples is simply to demonstrate how this approach works. Checks were made to quantitatively estimate sampling error and step-size error in τ . Substantial improvements in efficiency appear possible. It should be feasible to combine Ritz-Galerkin [84], spectral [77], or other methods [85] with boundary iteration. The MALA technique for accelerated convergence is also applicable [78].

F. Propagation with known end-points
The techniques described above can be used to calculate probabilities of a given path amongst all the possible quantum paths in phase-space. This requires the solution of a higher-dimensional PSDE. However, it is assumed a priori, that one already knows both the initial marginal x-distribution, P x (x 0 ), and the final marginal y distribution, P y (y f ). How can one proceed when it is necessary to know the Q-function forρ f at a future time?
While there are many possible approaches, here we outline some ways to achieve this within the stochastic framework: Ground states: To obtain a ground state or stationary state of finite entropy forĤ =Ĥ 1 +Ĥ 2 , one may proceed by adiabatic passage as in some experiments [86]. A stateρ 0 which is stationary forĤ 1 is constructed. This could be the non-interacting ground state. The full Hamiltonian is defined aŝ Here λ(t) is varied so that λ(0) = λ(2T ) = 0, and λ(T ) = 1. In the limit of slow passage, the dynamical path has known endpointsρ f =ρ 0 . The state at t = T is approximately stationary.
Transitional paths: If both initial and final distributions are known, samples of all intermediate paths for t 0 ≤ t ≤ t f can be calculated. This provides a means to understand quantum dynamical processes and the paradoxes of measurement theory, via the probability distribution of the trajectories that are sampled while reaching a known final quantum state. This is relevant to quantum ontology [10].
Dynamical solutions: To obtain a true dynamical solution, a known stateρ 0 at t 0 must be evolved to an unknown stateρ f at a time t f . This requires Metropolis or similar Monte-Carlo sampling, by viewing the dynamical equations as a means to generate samples of Q(φ 0 ). The algorithm involves constructing an initial estimated y f , together with an appropriate stochastic process for y f . Details of this procedure will be treated elsewhere.
Canonical ensembles: If a many-body state is known to be in a canonical ensemble at thermal equilibrium, then it is generally assumed thatρ = exp −β Ĥ − µN where β = 1/k B T , T is the temperature, and µ is the chemical potential. This can be handled through an 'imaginary time' calculation, such that dρ/dβ = − 1 2 Ĥ ,ρ + , which involves an anti-commutator rather than a commutator. The operator equation can be turned into a phase-space equation and treated in a similar to the dynamical case, with additional potential terms [87].
Transitional ensembles: If a canonical ensemble is known at two different values of both β and µ, then the stochastic techniques defined above can be used to define a transition path and evaluate transitional ensemble properties at other β and µ values.

Conditional measurements:
One or more future outputs y f may be the macroscopic result from an amplifier. If one measures this in the future, then the dynamics can be conditioned on knowing the value of y f , which may be used to infer information about an unknown state in the past.

V. EXAMPLES
Hamiltonians in quantum field theory of the type analyzed here usually have quadratic and quartic terms. In this section we consider several examples, with details in single-mode cases. Let the general Hamiltonian have the form H = H 0 + H S + H I . Here H 0 is a free field term, H S describes quadrature squeezing, found in Hawking radiation or parametric down-conversion, and H I is a quartic nonlinear particle scattering interaction.
Each of these cases will be treated separately below for simplicity, but they can be combined if required.

A. Free-field case
After discretizing on a momentum lattice, and using the Einstein summation convention, the free-field Hamiltonian can be written in normally-ordered form as The corresponding Q-function equations are: Hence, the coherent amplitude evolution equations are: The simplest case is a single-mode simple harmonic oscillator Hamiltonian, such that: H = ωâ †â . This corresponds to a characteristic equation ofα = −iωα. The expectation value of the coherent amplitude in the Qfunction has the equation: which is identical to the corresponding Heisenberg equation expectation value. There is no diffusive behavior or noise for these terms, and as a result the Q-function has an exactly soluble, deterministic quantum dynamics. The evolution is noise-free, with no need to make the transformations outlined above, since from (4.14), the steady-state in extra dimensions is given by solving(5.3). There is no difference here between classical and quantum dynamics, as pointed out by Schrödinger [88].

B. Squeezed state evolution
Next, we consider quadratic interaction terms that are mapped to second-order derivatives in the Q-function. These cause squeezed state generation and include quantum noise. They contain dynamics that leads to a model for quantum measurement as well as quantum paradoxes, including EPR and Bell inequality violations.
Following the notation of Eq (2.19), the general squeezing interaction term is H Such quadrature squeezing interactions are found in many areas of physics [89]. They illustrate how the Q-function equation behaves in the simplest nontrivial case where there is a diffusion term that is not positive-definite. We will investigate this in some detail, with numerical examples. This case illustrates very clearly how complementary variance changes are related to complementary time propagation directions.
Physically, these terms arise from parametric interactions, and lead to the dynamics that cause quantum entanglement. They are widespread, occurring in systems ranging from quantum optics to black holes, via Hawking radiation. The simplest case is a single-mode quantum squeezed state with

Q-function dynamics
We can calculate directly how the Q-function evolves in time. Applying the correspondence rules as previously, one obtains a Fokker-Planck type equation, now with second-order terms. Combining these terms into one equation gives: Hence, on using the real quadrature definitions of Eq (2.32) with e iη = i, and making a variable change so that iα = (x + iy) /2, we obtain This demonstrates the typical behavior of unitary Qfunction equations. The diffusion matrix it is traceless and equally divided into positive and negative definite parts. In this case the X + quadrature decays, but has positive diffusion, while the the X − quadrature shows growth and amplification, but has negative diffusion in the forward time direction. The amplified quadrature, which corresponds to the measured signal of a parametric amplifier, has a negative diffusion and therefore is constrained by a future time boundary condition.
If initially factorizable, the Q-function solutions can always be factorized as a product with Q = Q + Q − . Then, if t − = t 1 + t 2 − t, the time-evolution is diffusive, with an identical structure in each of two different time directions: The corresponding forward-backwards SDE is uncoupled, with decay and stochastic noise occurring in each time direction: where dw µ dw ν = 2δ µν dt. From these equations one can calculate immediately that: This equation for the variance time-evolution implies that the variance is therefore reduced in each quadrature's intrinsic diffusion direction, for an initial vacuum state, with the solution in forward time given by: x 2 (t) = 1 + e −2t y 2 (t) = 1 + e 2t . (5.11) Therefore, the variance reduction occurs in the forward time direction for x, giving rise to quadrature squeezing, and in the backward time direction for y, corresponding to gain in the forward time direction. However, neither anti-normally ordered variance is reduced below one. This is the minimum possible, corresponding to zero variance in the unordered operator case. With this choice of units, the diffusion coefficient is d = 2, so the overall Lagrangian is: The net effect of the stochastic processes in opposite time directions is that growth in the uncertainty of one quadrature in one time direction is cancelled by the reduction in uncertainty of the other quadrature in the opposite time direction. This behavior is shown in Figs (5) to (8), which illustrate numerical solutions of the forward-backward equations using the techniques of the previous section.
These solutions use 1600 trajectories, and hence include sampling error. Three dimensional graphs show equilibration in the extra dimension. Two dimensional graphs show results near equilibrium at τ = 5, with plots of variance in X ± vs time.

Comparison to operator equations
Defining quadrature operatorsŶ =â +â † andX = i â −â † , this physical system has the well-known behavior that variances change exponentially in time [90], in a complementary way. Given an initial vacuum state in which X 2 (0) = Ŷ 2 (0) = 1, the Heisenberg equa- Hence, theX quadrature is squeezed, developing a variance below the vacuum fluctuation level, and thê Y quadrature is unsqueezed, developing a large variance. This maintains the Heisenberg uncertainty product, which is invariant. Once operator ordering is taken into account, this gives an identical solution to the Q-function solution in Eq 5.11, because the operator correspondences are for antinormal ordering. If we use {} to denote this, then: In both cases there is a reduction in variance in the direction of positive diffusion. If there is an initial vacuum state, then quadrature squeezing occurs in X in the forward time direction, with a variance reduced below the vacuum level. Backward time squeezing occurs in Y , which also has forward-time gain.

Higher dimensional stochastic equation
In the matrix notation used elsewhere, this means that we have d = 2, and: with c = 0, so that the quantum dynamics occurs as the steady-state of the higher dimensional equation: with boundary values such that: These boundary conditions are known as mixed boundary conditions. They are partly Dirichlet (specified value), and partly Robin (specified linear combination of value and derivative). Numerical solutions for the the squeezed x equations are given in Figs (7) and (8), while those for the unsqueezed y equations are given in Figs (5) and (6).The effects of sampling error are seen through the two solid lines, giving one standard deviation variations from the mean. Exact results are included via the dashed lines.  Fig (4).

C. Quartic Hamiltonian example
While general quantum field Hamiltonians are certainly possible, here we treat the single-mode case to clearly illustrate the form of the relevant diffusion equation. This includes the most significant issues.
Following the notation of Eq (2.19), the single-mode nonlinear interaction term is: This single-mode problem can be solved using other methods, making it a useful benchmark [4]. However, when generalized to a nonlinear scalar field theory by including multiple modes and linear couplings, these analytic methods are no longer applicable. The addition of extra modes and linear couplings does not significantly change the arguments used here. We study this case in order to understand the effect of cross-couplings between the forward and backward time-evolution.

Fokker-Planck equation in complex phase-space
From the Q-function identities of (2.16) , after reordering the differential operators, and taking g = 1 for simplicity, one obtains: This demonstrates how the ordering identities apply. Damping and detuning terms are not included. For the quartic Hamiltonian, zeroth order and fourth order derivative terms cancel. This equation is known from earlier work in quantum optics [91]. As a simple check, one can integrate the Fokker-Planck equation in phase-space to obtain moments, hence showing that: Since the Q-function averages correspond to anti-normal ordering, one recovers the same expectation value dynamics as for the Heisenberg equations, which are:

Fokker-Planck equation with constant diffusion
We introduce a change of variable to a complex phase θ, with a scaling factor of i/2 to simplify the resulting algebra, so that: The result of changing variables in the distribution, is that in θ coordinates the distribution is modified by the Jacobian of the transformation, so that: One also must take account of the chain rule for derivatives when changing variables, which means that To transform to phase coordinates with constant diffusion, the Fokker-Planck equation is first transformed into a form that includes the effects of the Jacobian, followed by a variable change to the new variables. The combined effect this is that the equation for Q, after the variable change, is given by: where we have defined a number variable equivalent to the particle number as in (2.11), so that n ≡ αα * − 1.

D. Transformation to real coordinates
As proved in previous sections, in this equation the diffusion term is not positive definite. Accordingly, just as with the squeezing Hamiltonian, there is no equivalent forward time stochastic process. To show this, let θ = x + iy , and n = αα * − 1 = exp (x − y) − 1, so that: (5.26) As expected from the traceless diffusion property, the equation has a simultaneous positive diffusion in one real coordinate, and negative diffusion in the other, giving the result that: This means that the drift term is 28) and the forward and backwards equations are not factorizable, owing to the coupling term n (x, y), which is proportional to particle number. The total Lagrangian is: These equations are equivalent to a forward-backwards stochastic equation. The two stochastic equations are almost identical in each time direction, although with opposite drift terms: n(x (t ) , y(t ))dt + t t0 dw x y(t) = y(t f ) + Unlike the previous example, the two time directions are coupled to each other, since n depends on both fields. This implies that scattering takes place between the positive and negative time direction propagating fields. To solve for the quantum dynamical time evolution requires an understanding of the coupled evolution of both quadrature fields.
To obtain a dynamical solution from the coupled, forward-backward stochastic equations, we must transform this equation using the real, time-symmetric action principle. In this case, the equivalent extra-dimensional equation is: where ζ µ (t, τ ) ζ ν (t , τ ) = 2δ µν δ (τ − τ ) δ (t − t ). Thus, these extra-dimensional dynamical equations have a remarkably simple mathematical structure.

VI. SUMMARY
The existence of a time-symmetric probabilistic action principle for quantum fields has several ramifications. It describes a different approach to the computation of quantum dynamics. Neither imaginary time nor oscillatory path integrals are employed. More generally, time evolution through a symmetric stochastic action can be viewed as a dynamical principle in its own right. It is equivalent to the traditional action principle of quantum field theory. The advantage is that it is completely probabilistic, even for real-time quantum dynamics.
An interesting property of this method is that it can provide an ontological interpretation of quantum mechanics. This is described in greater detail elsewhere. The picture is that physical fields can propagate both from the past and from the future. This is different to a hidden variable theory. These only allow causality from past to future. As a result, one can have quantum features including vacuum fluctuations, sharp eigenvalues and Bell violations [10].
The power of rapidly developing petascale and exascale computers appears well-suited to these approaches. Enlarged spatial lattices and increased parallelism are certainly needed. Yet this may not be as problematic to handle as either exponential complexity or the phase problems that arise in other approaches. It is intriguing that the utility of an extra dimension is widely recognized both in general relativity and quantum field theory. One may speculate that extending this action principle to curved space-time may yield novel quantum theories. This could lead to new approaches to unification.

ACKNOWLEDGMENTS
PDD thanks the hospitality of the Institute for Atomic and Molecular Physics (ITAMP) at Harvard University, supported by the NSF, and the Weizmann Institute of Science through a Weston Visiting Professorship. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. It was also funded through an Australian Research Council Discovery Project Grant DP190101480.

APPENDIX
In this Appendix the higher-dimensional equilibration equations are obtained for a more general case, with a diffusion term d µ (t) that is time and index dependent. We recall that G ([φ], τ ) satisfies a higher-dimensional functional partial differential equation of: To obtain the steady-state solution as τ → ∞, (3.32) for G, from Eq (3.33) A[φ] must satisfy: L c (φ,φ, t)dt , (6.2) where the central-difference Lagrangian is: This leads to Euler-Lagrange type equations with a drift defined as: where Solving for the higher-dimensional drift, the timederivative term is: and the remaining potential term is: Introducing first and second derivatives,φ ≡ ∂φ/∂t andφ ≡ ∂ 2 φ/∂t 2 , there is an expansion for the higherdimensional drift term A in terms of the field timederivatives: A =φ + cφ + a . (6.8) Here, c is a circulation matrix, while a is a pure drift: Here partial derivatives with respect to time indicate derivatives for the explicitly time-dependent terms only, where the Hamiltonian coefficients are changing in time. The function U is an effective potential, which acts to generate an effective force on the trajectories: The functional Fokker-Planck equation given above is then equivalent to a stochastic partial differential equation (SPDE): where the stochastic term ζ is a real deltacorrelated Gaussian noise such that ζ µ (t, τ ) ζ µ (t , τ ) = 2d µ (t) δ µν δ (t − t ) δ (τ − τ ). The final stochastic partial differential equation in τ that φ must satisfy in detail is then: ∂φ ∂τ =φ + cφ + a + ζ (t, τ ) . (6.12)