Quantum mechanics and neural networks

When discussing a Euclidean field theory in d spacetime dimensions, the local observables of fundamental interest are the correlation functions

$$\begin{align} G^{\left(n\right)} \left( y_1, \ldots, y_n \right) = \langle \phi \left( y_1 \right) \ldots \phi \left( y_n \right) \rangle \, ,\end{align} \tag{ 2.1 }$$

where φ is a local operator in the spectrum of the theory (which may be a scalar or carry additional indices that are suppressed) and the quantities $y_i \in \mathbb{R}^d$, for $i = 1, \ldots, n$, are spacetime points at which the field operators have been inserted. The correlators $G^{(n)}$ are also referred to as Schwinger functions. These are the primary objects of study in constructive approaches to QFT, where one investigates whether these functions satisfy certain conditions such as the OS axioms, which guarantee that a Euclidean theory may be analytically continued to a corresponding Lorentzian theory. When $d \unicode{x2A7E} 2$, it is natural to refer to the theory associated with such correlation functions as a field theory, since the degree of freedom φ varies in both Euclidean time and in space.

Restricting to the case d = 1 that will be relevant for quantum mechanics, we instead consider ‘fields’ — such as the position x(t) of a quantum particle subject to a harmonic potential—which depend on a time variable $t \in \mathbb{R}$, but not on any spatial coordinates. Therefore, in a quantum mechanical theory, the natural observables are instead correlators

$$\begin{align} G^{\left(n\right)} \left( t_1 , \ldots , t_n \right) = \langle x \left( t_1 \right) \ldots x \left( t_n \right) \rangle \, ,\end{align} \tag{ 2.2 }$$

which depend on n time coordinates $t_i \in \mathbb{R}$. We view the correlation function (2.2) as the expectation value, in some probability distribution, of a product of random variables $x ( t_i )$ for $i = 1 , \ldots , n$. A theory of quantum mechanics should yield a prescription for computing any such expectation value, for any choice of times t_i and for any finite n.

However, not all collections of correlation functions $G^{(n)}$ are suitable to be interpreted as arising from physically realistic quantum theories. For instance, we envision the random variables $x ( t_i )$ as representing points on the trajectories of quantum particles. On physical grounds, any such trajectory should be a continuous function of Euclidean time, since a particle should not be able to ‘jump’ instantaneously from one position to another. We will therefore impose a technical condition called ‘mean-square continuity’ on the distribution of random variables $x(t_i)$ in our quantum model. Formally, we require that

$$\begin{align} \lim_{t \to s} \left\langle \, \left| x \left( t \right) - x \left( s \right) \right|^2 \, \right\rangle = 0\end{align} \tag{ 2.3 }$$

for all s. In particular, this condition implies that the sample paths x(t) are continuous and that the two-point function $G^{(2)} ( t, s )$ is a continuous function of t and s⁹.

A second basic requirement comes from the Kallén–Lehmann spectral representation [24, 25] of the two-point function, which is a non-perturbative result that holds for quite general quantum field theories. In one Euclidean spacetime dimension, this representation for the two-point function $G^{(2)} ( t, s )$ in any theory of quantum mechanics takes the form

$$\begin{align} G^{\left(2\right)} \left( t, s \right) = \int_0^{\infty}\mathrm d m \, \rho \left( m \right) \mathrm e^{-m | t - s |} \, ,\end{align} \tag{ 2.4 }$$

where $\rho ( m )$ is a finite Borel measure, which means that ρ is positive definite and that its integral is finite (though not necessarily equal to 1, in which case it is called a probability measure). To see why this assumption is reasonable, it is convenient to consider a conventional quantum system with a discrete spectrum that is bounded below. After shifting the Hamiltonian by a constant so that the ground state energy is zero, one may insert a complete set of energy eigenstates $1 = \sum_n \ket{n} \bra{n}$ to find

$$\begin{align} \langle x \left( t \right) x \left( s \right) \rangle = \sum_{n = 1}^{\infty} \mathrm e^{- E_n | t - s | } \left| \langle 0 \mid x \mid n \rangle \right|^2 = \int_0^{\infty}\mathrm d m \, \rho \left( m \right) \mathrm e^{-m | t - s |} \, ,\end{align} \tag{ 2.5 }$$

where

$$\begin{align} \rho \left( m \right) = \sum_{n = 1}^{\infty} \left| \langle 0 \mid x \mid n \rangle \right|^2 \delta \left( m - E_n \right) \, .\end{align} \tag{ 2.6 }$$

Although this result is most often applied in higher-dimensional field theory, the version of this expansion for 1d Euclidean systems was already used in some early works such as [26]. Since $\rho ( m )$ must be normalizable, equation (2.4) implies that the two-point function $G^{(2)} ( t, s )$ must be finite for any values of t and s. Indeed, in the coincident-point limit one has $G^{(2)} ( t, t ) = \int_0^\infty\mathrm d m \, \rho ( m ) \lt \infty$, and for any other t ≠ s the integrand is positive-definite and upper-bounded by $\rho ( m )$, so $G^{(2)} ( t, s )$ converges to a finite value for any $t, s$.

These two basic physical inputs are encoded in the following two properties.

• (P1)  
  The expectation values in our quantum model obey the mean-square continuity condition (2.3). This implies that the two-point¹⁰ function $G^{(2)} ( t, s )$ is a continuous function of its arguments, and that the sample paths x(t) are continuous in t.
• (P2)  
  The two-point function $G^{(2)} ( t, s )$ is finite for all values of the input times t, s.

We therefore adopt the definition that a ‘minimal quantum mechanical model’ or MQM is an assignment, for any $n \in \mathbb{N}$ and to every collection of points $t_1 , \ldots, t_n$, of a joint probability distribution for a collection of random variables $x ( t_1 )$, $\ldots$, $x ( t_n )$, such that the expectation values in this distribution obeys (P1) and (P2).

Here we have focused on Schwinger functions $G^{(n)}$ rather than other common quantities of interest in quantum theories, such as the eigenvalues of the Hamiltonian. However, given the collection of all such correlation functions (2.2) for arbitrary times, it is possible to reconstruct other desired observables. For instance, the energy spectrum can be extracted by performing exponential fits of the fall-offs of correlation functions at large time separations. In fact, for quantum systems that satisfy the OS axioms (which we will discuss in section 3) and can thus can be continued to real time, the OS reconstruction theorem guarantees that the correlation functions $G^{(n)}$ are sufficient to reconstruct the Hilbert space and all local observables of the real-time theory [14]. Therefore, at least for models obeying the OS axioms, we suffer essentially no loss of generality by restricting attention to correlation functions.

Furthermore, although here we primarily discuss quantum systems with a single degree of freedom, one can straightforwardly generalize this definition to models with a collection of bosonic coordinates $x^i ( t )$, to models with anticommuting fields $\psi^i ( t )$ (such as supersymmetric quantum mechanics), and to other settings. For instance, this allows us to represent many-body quantum systems, where the index on the coordinates $x^i ( t )$ label degrees of freedom associated with distinct sites or particles. Let us also remark that seemingly ‘discrete’ quantum models, such as a two-level system, can be accommodated in this framework. For instance, to describe a qubit, we choose three random variables $x^i ( t ) = \left( \phi ( t ) , \theta ( t ) , \psi ( t ) \right)$ using the Euler angle representation, as in the standard construction of the path integral for a spin-$\frac{1}{2}$ particle (see e.g. section 3.3 of [27] for a review). In these more general settings, one could engineer additional symmetries beyond those considered in this work, such as invariance under rotations $x^i_t \to {R}{^{i}_{j}} x_{j, t}$ or supersymmetry transformations mixing bosons and fermions; one approach to doing this is to use the ‘parameter absorption trick’ described in [28].

We refer to the models defined above as ‘minimal’ because, depending on which properties one insists upon as necessary for a theory to qualify as genuinely quantum, these conditions may not be sufficient, although we believe it is fair to say that they are necessary. For example, one might reasonably adopt the view that a Euclidean (imaginary time) model should only be called ‘quantum’ if it can be analytically continued to a Lorentzian (real time) model with unitary time evolution. We will eventually demand that our systems satisfy additional properties, such as time translation symmetry and RP, which ensure the existence of such a Lorentzian continuation.

Our definition of a minimal quantum system, obeying (P1) and (P2), is a special case of the mathematical notion of a SP. In modern treatments, the general definition of a SP is phrased in terms of a σ-algebra $\mathcal{F}$ of events in a set Ω called the sample space, along with a measure $\mathbb{P}$ which assigns probability to elements of $\mathcal{F}$. However, for our purposes it will suffice to work with a more elementary definition: by a SP, we mean a collection of real-valued random variables

$$\begin{align} \left\{x \left( t \right) \mid t \in T \right\} \, ,\end{align} \tag{ 2.7 }$$

where T is the index set for the process, which we typically choose to be either $\mathbb{R}$ or a finite interval $[a, b]$. We will often abbreviate a SP x(t) with a subscript, as in x_t. Since for each fixed $t_i \in \mathbb{R}$ the quantity $x ( t_i )$ is a random variable, one can speak of the joint probability distribution¹¹

$$\begin{align} \mathbb{P} \left( x \left( t_1 \right) , \ldots , x \left( t_n \right) \right)\end{align} \tag{ 2.8 }$$

for any finite collection of these random variables. The set of all such joint probability distributions, for any collection of n variables t_i in the index set T and for any n, is the defining data of the SP x_t.

We therefore see that the definition of a minimal quantum theory given above is an example of a SP, since the defining data for both of these objects is an assignment of a joint probability distribution (2.8) to every finite collection of times t_i. The only distinction between these notions is that a minimal quantum theory must satisfy conditions (P1) and (P2), which means that such minimal quantum models are a proper subset of SPs. More precisely, a MQM is a SP x_t which is mean-square continuous, i.e. it obeys (2.3), and square-integrable, which implies that the correlation function $G^{(2)} ( t, s )$ is finite for every pair of inputs t and s, ensuring (P2).

The preceding definition referred only to the values $x ( t_i )$ of the process x_t at specific points $t_i \in \mathbb{R}$. However, a complementary viewpoint of a SP is as a probability distribution over the space of functions $x : \mathbb{R} \to \mathbb{R}$. From this perspective, one performs a random ‘draw’ of a function x_t, and then the values of this sample function at points t_i furnish us with realizations of the random variables $x ( t_i )$. This interpretation of a SP is aligned with our usual intuition in the path integral formulation of quantum mechanics, where we view expressions like

$$\begin{align} \langle x \left( t_1 \right) x \left( t_2 \right) \rangle = \int \mathcal{D} x \, x \left( t_1 \right) x \left( t_2 \right) \mathrm e^{- S \left[ x \right]}\end{align} \tag{ 2.9 }$$

as formally computing an expectation value with respect to a probability measure on the space of paths which is determined by the action $S[x]$ and path integral measure $\mathcal{D} x$. However, we note that the SPes (or minimal quantum systems) that we consider here are not assumed to be defined by an action: instead we work directly with an abstract probability distribution (2.8) which plays the role of the combination $\mathcal{D} x \, \mathrm e^{- S [ x ] }$ and thus captures the combined information (both classical and quantum, whereas the action $S[x]$ alone would define only the classical dynamics) of the model.

In some mathematical treatments of SPs, one avoids this picture of a ‘random function’ (which requires imposing additional regularity assumptions on x_t for technical reasons), preferring to speak only of the random variables $x ( t_i )$ associated with specific points in time. We will not concern ourselves with such technical subtleties in this work, and instead use the two perspectives interchangeably.

If one wishes to take this interpretation of a SP x_t as a random function seriously, it is natural to ask whether one can parameterize the family of functions that are randomly generated by the process x_t in terms of some ‘coordinates’ which are real-valued quantities that label a particular function in the family. This leads us to the other main definition of interest in this work, which is that of a NN.

For the purposes of this work¹², a NN is a parameterized family of functions $\phi_\theta : \mathbb{R} \to \mathbb{R}$ determined by a set of parameters θ, which are random variables with a joint probability distribution $P ( \theta )$. There can, in general, be many such parameters θⁱ for $i = 1 , \ldots, N$ (or infinitely many), but we will suppress indices and write θ to indicate the collection of all parameters. Performing a random draw of the parameters θ then gives a specific instance of a random function φ_θ. A simple example of a NN is

$$\begin{align} \phi_\theta \left( t \right) = \cos \left( t + \theta \right) \, , \quad \theta \sim U \left[ 0 , 2 \pi \right] \, ,\end{align} \tag{ 2.10 }$$

where the latter notation indicates that θ is drawn from a uniform distribution on the interval $[0 , 2 \pi]$. We see that the function φ_θ therefore inherits randomness from the random variable θ; a particular draw of the parameter, say $\theta = \pi$, then determines an instance of the NN, here $\phi_\pi ( t ) = \cos ( t + \pi )$. We refer to the fixed functional form that describes the dependence of φ_θ on its parameters, such as the cosine function of (2.10), as the architecture of the NN, and we refer to the probability distribution $P ( \theta )$ over the parameters θ as the parameter density.

Clearly every NN satisfies the definition of a SP, since taken together, the architecture and joint probability distribution $P ( \theta )$ over the parameters determine a joint probability distribution

$$\begin{align} \mathbb{P} \left( \phi_\theta \left( t_1 \right) , \ldots, \phi_\theta \left( t_n \right) \right)\end{align} \tag{ 2.11 }$$

for the collection of random variables $\phi_\theta ( t_i )$ given any finite collection of times t_i. We say that a SP x_t which admits a representation using a NN, so that $x_t = \phi_\theta ( t )$ with some parameter density $P ( \theta )$, is a NN SP or NN-SP; similarly, any minimal quantum model $x_t = \phi_\theta ( t )$ which can be represented as a NN is said to be a NN quantum mechanics or NN-QM. Said differently, a NN-QM is a NN-SP that enjoys properties (P1) and (P2). For any NN-SP, the correlation functions (2.2) can be computed using parameter space integrals as

$$\begin{align} G^{\left(n\right)} \left( t_1 , \ldots, t_n \right) = \int\mathrm d \theta \, P \left( \theta \right) \, \phi_\theta \left( t_1 \right) \ldots \phi_\theta \left( t_n \right) \, ,\end{align} \tag{ 2.12 }$$

where the integral runs over all of the parameters θ but again we suppress indices.

Given that every NN $\left( \phi_\theta , P ( \theta ) \right)$ determines a SP, one might ask about the opposite direction: does every SP x_t admit a NN representation in terms of an architecture φ_θ and a parameter density $P ( \theta )$? One could answer this question in a rather trivial way by choosing an uncountably infinite set of parameters θ^t, indexed by $t \in \mathbb{R}$, and then declaring that $\phi_\theta ( t ) = \theta^t$. This is simply choosing a separate random variable for the output of the NN at each point. Since the data of the SP x_t involves a joint probability distribution (2.8) for each collection of points t_i, one could then define the parameter density $P ( \theta )$ to agree with the joint distribution of x_t on any finite set of inputs. However, this construction is unwieldy and not very interesting.

One might instead ask whether every SP can be represented via a NN in a useful way, with only countably many parameters and in terms of an architecture which is amenable to theoretical analysis. The answer to this question is provided by a theorem discussed by Kosambi in [29], by Karhunen in [30], and by Loève in [31], which we now recall.

Theorem 2.1. (Kosambi–Karhunen–Loève). Let x_t be a square-integrable SP defined on an interval $[a, b]$ such that $\langle x ( t ) \rangle = 0$ for all $t \in [ a, b ]$, and suppose that the two-point function $\langle x ( t_1 ) x ( t_2 ) \rangle$ is continuous. Then x_t admits a decomposition

$$\begin{align} x_t = \sum_{k = 1}^{\infty} \theta^k e_k \left( t \right) \, ,\end{align} \tag{ 2.13 }$$

where e_k is a set of continuous, orthogonal real-valued functions on $[a, b]$ and θ^k are a collection of pairwise uncorrelated (but not necessarily¹³ independent) random variables.

First let us discuss the assumptions of the KKL theorem. Two of the conditions of the theorem—continuity of the two-point function and square-integrability—are implied by (P1) and (P2), which are satisfied by any MQM by definition. The last requirement of theorem 1 is that the SP x_t has zero expectation value for all t, which is not a condition which we impose in general. However, this is easy to remedy. Consider some x_t which is square-integrable and has a continuous two-point function but which does not satisfy the zero-mean condition. Any square-integrable SP has finite first moment $\langle x ( t ) \rangle$ for all t, by the Cauchy–Schwarz inequality. We define

$$\begin{align} \hat{x}_t = x_t - \langle x \left( t \right) \rangle \, .\end{align} \tag{ 2.14 }$$

Then $\hat{x}_t$ satisfies all of the assumptions of the Kosambi–Karhunen–Loève theorem. It follows that the original SP can be written as the sum of the KKL decomposition for $\hat{x}_t$ and the deterministic function $\langle x ( t ) \rangle$. Viewing the expansion coefficients θ^k appearing in (2.13) as NN parameters, and choosing the architecture

$$\begin{align} \phi_\theta \left( t \right) = \langle x \left( t \right) \rangle + \sum_{k = 1}^{\infty} \theta^k e_k \left( t \right) \, ,\end{align} \tag{ 2.15 }$$

we see that every square-integrable SP x_t with continuous two-point function admits a NN representation, albeit generically one with a countably infinite set of parameters. As we mentioned, every minimal quantum model obeys (P1) and (P2) and thus satisfies the conditions of square-integrability and continuity of $G^{(2)}$, so we conclude that every MQM is also a NN-QM. In fact, it is now redundant to use the term ‘NN-QM’ since we have proven that the space of NN-QMs (i.e. SPs obeying (P1) and (P2)) precisely coincides with the space of MQMs. A Venn diagram of various relevant subsets of the set of SPs appears in figure 1.

Zoom In Zoom Out Reset image size

This observation furnishes us with a sort of ‘quantum universal approximation theorem,’ in the sense that any MQM can be approximated to arbitrary accuracy by a truncated architecture of the form (2.15) which includes finitely many of the terms in the sum¹⁴. As the number of included terms grows larger, this architecture becomes a better and better approximation to the true statistics of the SP. The specific expansion generated by the Kosambi–Karhunen–Loève theorem, where the $e_k ( t )$ are chosen to be eigenfunctions of a certain linear operator associated with the two-point function of the process, is ‘optimal’ in the sense that a truncated KKL expansion gives the smallest mean-square error among all truncated orthonormal expansions of the SP. The result (2.13) is primarily a representation theorem, analogous to universal approximation theorems for continuous functions, and in practice one might prefer to represent a MQM in a different way which is more amenable to training, such as using a deep NN. We will return to this point in section 4.

Let us emphasize that this result holds for any quantum model obeying (P1) and (P2), and is not restricted to cases which satisfy additional conditions that we typically assume when modeling most physical quantum systems, such as time translation invariance or RP (again, the latter leads to unitary time evolution). This is important since each of these standard assumptions about quantum systems is relaxed in some use cases. For instance, in the study of Floquet systems, one does not impose continuous time translation symmetry, instead assuming the weaker condition that the Hamiltonian is periodic in time, i.e. $H ( t + T ) = H ( t )$ for some period T. Likewise, when considering open quantum systems, one does not generally have unitary time evolution; rather the dynamics are described by the Lindblad master equation, which models a system that interacts with some environment or ‘bath’ and thus evolves in a non-unitary way. Both of these applications can be accommodated within the framework of NN-QM, since we have made only the mild assumptions (P1), (P2) about the underlying process x_t.

Having said this, we will now set aside these ‘unconventional’ applications, and for the remainder of this work we will focus on standard quantum systems which enjoy time translation invariance and unitary time evolution. The additional assumptions on x_t, or φ_θ, which guarantee these properties are the subject of the next section.

In this section, we will pass to the NN description of a quantum system and outline one mechanism by which RP (E4) can be achieved. The advantage of the NN formulation is that we may write expectation values as explicit parameter space integrals, which can make certain manipulations easier.

Consider a NN quantum system x_t which is presented via an architecture $x ( t ) = \phi_\theta ( t )$ with associated parameter density $P ( \theta )$. First we will express the condition for RP in this language. Fix a collection of positive times $t_1 , \ldots, t_n \gt 0$ and a (possibly complex) function F of n variables. The condition (3.1) then reads

$$\begin{align} \int \mathrm d \theta \, P \left( \theta \right) \, F \left( \phi_\theta \left( t_1 \right) , \ldots, \phi_\theta \left( t_n \right) \right) \left( F \left( \phi_\theta \left( - t_1 \right) , \ldots, \phi_\theta \left( - t_n \right) \right) \right)^\ast \unicode{x2A7E} 0 \, ,\end{align} \tag{ 3.4 }$$

where again $\ast$ denotes complex conjugation. RP requires that the inequality (3.4) holds for any suitable function F and any finite collection of positive times.

Let us now ask how one might engineer a NN-QM such that RP is automatically satisfied. One obvious way to ensure that a quantity is positive is to write it as a perfect square. We can attempt to force the integrand of (3.4) to be a perfect square by ‘splitting’ the parameter space integral and writing it in a factorized form as follows. Assume that the set $\left\{\theta \right\}$ of parameters associated with the NN-QM $( \phi_\theta , P ( \theta ) )$ decomposes into three mutually disjoint subsets:

$$\begin{align} \left\{\theta \right\} = \left\{\theta^0 \right\} \cup \left\{\theta^+ \right\} \cup \left\{\theta^- \right\} \, .\end{align} \tag{ 3.5 }$$

As usual, there may be many parameters in each subset, but we suppress indices and write θ⁰, $\theta^+$, and θ⁻ to denote the collection of all parameters in each subset.

The critical assumption we now make is the following. Suppose that, for positive times t > 0, the architecture $\phi_\theta ( t )$ depends on the parameters $\theta^+$ and θ⁰ but not on θ⁻. Likewise, at negative times t < 0, assume that $\phi_\theta ( t )$ depends on θ⁻ and θ⁰ but not $\theta^+$. This condition allows us to rewrite (3.4) as

$$\begin{align} \langle F \left( T F \right)^\ast \rangle & = \int\mathrm d \theta^0 \, P \left( \theta^0, \theta^+, \theta^- \right) \, \left( \int\mathrm d \theta^+ F \left( \phi_\theta \left( t_1 \right) , \ldots, \phi_\theta \left( t_n \right) \right) \right) \nonumber \\ &\quad \quad \cdot \left( \int\mathrm d \theta^- \left( F \left( \phi_\theta \left( - t_1 \right) , \ldots, \phi_\theta \left( - t_n \right) \right) \right) \right)^\ast \, ,\end{align} \tag{ 3.6 }$$

where we have introduced the shorthand $\langle F \left( T F \right)^\ast \rangle$ for the integral on the left side of (3.4). Here T is the time-reversal operator which acts on any function of a collection of insertions $\phi_\theta ( t_i )$ by reversing each time as $T : \phi_\theta ( t_i ) \mapsto \phi_\theta ( - t_i )$.

We have succeeded in factorizing the F and $\left( T F \right)^\ast$ insertions, but not the parameter density. To achieve the latter, let us pause to consider which conditions on $P ( \theta^0, \theta^+, \theta^- )$ are natural to impose. Although we are presently focused on RP, we are ultimately interested in engineering a NN-QM that obeys all of the OS axioms. In particular, if our NN satisfies the OS axiom (E1), then the associated SP $x_t = \phi_\theta ( t )$ is symmetric and hence x_t and x_−t are equivalent processes. This equivalence means that various joint probability distributions must be equal, such as

$$\begin{align} \mathbb{P} \left( \phi_\theta \left( t_1 \right) , \ldots , \phi_\theta \left( t_n \right) \right) = \mathbb{P} \left( \phi_\theta \left( - t_1 \right) , \ldots , \phi_\theta \left( - t_n \right) \right) \, .\end{align} \tag{ 3.7 }$$

In the NN description, joint probability distributions on outputs such as (3.7) are controlled by the parameter density. A sufficient set of conditions that one may impose on the density $P ( \theta )$ which ensures that $\phi_\theta ( t )$ generates a symmetric SP is

$$\begin{align} P \left( \theta^0, \theta^+, \theta^- \right) = P_0 \left( \theta^0 \right) \cdot P_+ \left( \theta^0, \theta^+ \right) \cdot P_- \left( \theta^0, \theta^- \right) \, ,\end{align} \tag{ 3.8 }$$

along with an assumption that time-reversal of the input t can be ‘absorbed’ into an action on the parameters, in the sense that for t > 0 the architecture obeys

$$\begin{align} \phi_{\theta^0 , \theta^+} \left( t \right) = \phi_{\theta^0 , \theta^{+ ^{\prime}}} \left( - t \right) \, , \quad \phi_{\theta^0 , \theta^-} \left( - t \right) = \phi_{\theta^0 , \theta^{- ^{\prime}}} \left( t \right) \, ,\end{align} \tag{ 3.9 }$$

where $\theta^{\pm ^{\prime}}$ are re-defined sets of parameters, such that

$$\begin{align} P_+ \left( \theta^0 , \theta^{- ^{\prime}} \right) = P_- \left( \theta^0 , \theta^- \right) \, , \quad P_- \left( \theta^0 , \theta^{+ ^{\prime}} \right) = P_+ \left( \theta^0 , \theta^+ \right) \, .\end{align} \tag{ 3.10 }$$

Equation (3.10) looks somewhat strange due to expressions like $P_\pm \left( \theta^0 , \theta^{\mp ^{\prime}} \right)$ which involve evaluation of $P_\pm$ on inputs with the ‘wrong-sign’ subscripts $\mp$. However, let us emphasize that all of the parameters in these equations are simply dummy variables. At the risk of being pedantic, suppose there are N₀ parameters θ⁰ and $N_{\pm}$ parameters $\theta^{\pm}$. We assume that $N_+ = N_- = N$, so that both $P_+$ and $P_-$ are functions of $N_0 + N$ variables. Then the first equation of (3.10) states that the output of the function $P_-$, when evaluated on the $N_0 + N$ inputs θ⁰ and θ⁻, matches the output of the (different) function $P_+$, when $P_+$ is evaluated on the (different) set of $N_0 + N$ inputs θ⁰ and $\theta^{- ^{\prime}}$. Here $\theta^{- ^{\prime}}$ is a new set of N variables that are obtained by applying some $\mathbb{Z}_2$ transformation to the N variables θ⁻.

This ‘parameter absorption mechanism’ was first introduced in [4] to engineer NN field theories with particular symmetries. The specific mechanism of equations (3.9)–(3.10) has been chosen so that the NN-QM defined by $x_t = \phi_\theta ( t )$ enjoys time-reflection invariance, i.e. so that x_t and x_−t are equivalent SPs.

Finally, let us return to RP. With the additional conditions on the parameter density described above, (3.6) becomes

$$\begin{align} \langle F \left( T F \right)^\ast \rangle & = \int\mathrm d \theta^0 \, P \left( \theta^0 \right) \, \left( \int\mathrm d \theta^+ \, P_+ \left( \theta^0 , \theta^+ \right) \, F \left( \phi_\theta \left( t_1 \right) , \ldots, \phi_\theta \left( t_n \right) \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm d \theta^- \, P_- \left( \theta^0 , \theta^- \right) \, \left( F \left( \phi_\theta \left( - t_1 \right) , \ldots, \phi_\theta \left( - t_n \right) \right) \right) \right)^\ast \, \nonumber \\ & = \int\mathrm d \theta^0 \, P \left( \theta^0 \right) \, \left( \int\mathrm d \theta^+ \, P_+ \left( \theta^0 , \theta^+ \right) \, F \left( \phi_\theta \left( t_1 \right) , \ldots, \phi_\theta \left( t_n \right) \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm d \theta^- \, P_- \left( \theta^0 , \theta^- \right) \, \left( F \left( \phi_{\theta^{^{\prime}}} \left( t_1 \right) , \ldots, \phi_{\theta^{^{\prime}}} \left( t_n \right) \right) \right) \right)^\ast \nonumber \\ & = \int\mathrm d \theta^0 \, P ( \theta^0 ) \, \left( \int\mathrm d \theta^+ \, P_+ ( \theta^0 , \theta^+ ) \, F \left( \phi_\theta ( t_1 ) , \ldots, \phi_\theta ( t_n ) \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm d \theta^{- ^{\prime}} \, P_+ ( \theta^0 , \theta^{- ^{\prime}} ) \, \left( F \left( \phi_{\theta^{^{\prime}}} ( t_1 ) , \ldots, \phi_{\theta^{^{\prime}}} ( t_n ) \right) \right) \right)^\ast \nonumber \\ & = \int\mathrm d \theta^0 \, P ( \theta^0 ) \, \left| \int\mathrm d \theta^+ \, P_+ ( \theta^0 , \theta^+ ) \, F \left( \phi_\theta ( t_1 ) , \ldots, \phi_\theta ( t_n ) \right) \right|^2 \, .\end{align} \tag{ 3.11 }$$

In the first step, we have used (3.9) to absorb all minus signs in the evaluations of φ_θ at negative times into redefinitions of the parameters. We then applied (3.10) to replace $P_-$ in the second integral with $P_+$ evaluated on the transformed parameters, and we have changed integration variables from θ⁻ to $\theta^{- ^{\prime}}$, which we assume gives a trivial Jacobian so that no additional factors are generated from the change of measure¹⁸. Finally, in the last step we recognize that both $\theta^+$ and $\theta^{- ^{\prime}}$ are dummy variables that are integrated over, which means that the integrals inside the two sets of parentheses are in fact identical. This allowed us to write $\langle F \left( T F \right)^\ast \rangle$ as an integral of a non-negative density $P ( \theta^0 )$ multiplied by a perfect square, which is thus also non-negative, as desired. We conclude that

$$\begin{align} \langle F \left( T F \right)^\ast \rangle \unicode{x2A7E} 0\end{align} \tag{ 3.12 }$$

for a NN-QM satisfying the assumptions described above. This establishes RP for any quantum system represented by a NN that enjoys parameter splitting.

Let us now summarize the (rather involved) argument presented above, and make some further comments. We have determined a set of sufficient conditions for a NN-QM to enjoy RP. The key ingredients in our mechanism are:

• (I)  
  The parameter-space degrees of freedom which are relevant for positive times and negative times decouple from one another. By this, we mean that there is one set of parameters $\theta^+$ that describe the NN-QM at positive times, but which have no effect on the model at negative times; likewise, there is another set of parameters θ⁻ that affect the NN at negative times but not positive times.
• (II)  
  The parameter density enjoys a factorized form (3.8), and one can employ a ‘parameter absorption trick’ to relate the densities at positive and negative times. This assumption also implies time-reflection symmetry (not to be confused with RP).

The first ingredient (I) is reminiscent of locality in time; the behavior of the NN at large positive times is independent of its behavior in the far past, at negative times, which suggests some form of local interactions.

Alternatively, one can view the parameter splitting as a sort of ‘parameter space Markov property’ in the following sense. Suppose we view the parameters $\theta^+$ as describing $\phi_\theta ( t )$ for $t \unicode{x2A7E} 0$, θ⁻ as describing the NN for $t \unicode{x2A7D} 0$, and θ⁰ as determining the behavior of $\phi_\theta ( 0 )$. Say that we are given the output of a fixed realization of the NN on the half-line $t \unicode{x2A7D} 0$ and wish to predict its outputs in the future, at positive times. Parameter splitting roughly states that, to make this prediction, it is not necessary to condition on the NN outputs at negative times (which depend on θ⁻), but only on the most recent time t = 0 associated with θ⁰. This can be interpreted as a ‘memorylessness’ assumption, at least around the point t = 0. The connection between Markov-like behavior and RP is not an accident; in section 3.2 we will return to this topic, although there we consider the more conventional ‘output-space Markov property’ rather than this ‘parameter-space Markov property’ represented by the splitting mechanism.

We also note that, while our parameter splitting mechanism is sufficient for RP, it depends on the precise parameterization of the NN. A generic quantum system can be presented by NNs in multiple equivalent ways, and even if one representation enjoys parameter splitting, another presentation may fail to split. Therefore it cannot be necessary that a NN with a given parameterization must enjoy a parameter splitting in order to be RP. A simple way to see this is to consider the KKL decomposition of a minimal quantum mechanical model on an interval $[-T, T]$, which we repeat:

$$\begin{align} \phi_\theta \left( t \right) = \langle x \left( t \right) \rangle + \sum_{k = 1}^{\infty} \theta^k e_k \left( t \right) \, .\end{align} \tag{ 3.13 }$$

Each of the continuous, pairwise orthogonal functions $e_k ( t )$ on $[-T, T]$ generically has support both on the positive half-interval $[0, T]$ and on the negative half-interval $[-T, 0]$. For instance, in the KKL decomposition of a Brownian path or an Ornstein–Uhlenbeck process, the $e_k ( t )$ can be chosen to be trigonometric functions. In such a parameterization, there is no splitting of the parameters θ^k into appropriate subsets $\theta^+$ and θ⁻; every one of the θ^k affects the outputs of the NN φ_θ at both positive and negative times. But on the other hand, every MQM admits a KKL representation, including the reflection-positive processes. Therefore an RP process, presented via a NN with some choice of parameterization, need not exhibit parameter splitting with respect to that particular parameterization. Said differently, the existence or non-existence of a parameter splitting is a question which depends on the given NN representation $\phi_\theta ( t )$ for a quantum system, rather than an intrinsic question about the system itself¹⁹.

Engineering parameter splitting about t = 0

The preceding subsection explained that a set of parameter splitting assumptions on a NN-QM implies RP. One might then ask whether this set of assumptions is actually satisfied in any explicit examples. As it turns out, given any NN $\left( \phi_\theta , P ( \theta ) \right)$, it is possibly to transform this NN into a modified model $\left( \widetilde{\phi}_\theta , \widetilde{P} ( \theta ) \right)$ which obeys the parameter splitting requirements and therefore exhibits RP. In principle, this allows us to generate a large class of examples of NN-QM systems which enjoy RP, although they generally do not satisfy the other OS axioms.

Let us now describe the construction. Begin with any NN-QM defined by an architecture $\phi_\theta : \mathbb{R} \to \mathbb{R}$ and a probability distribution $P ( \theta )$ over parameters. In addition to the parameters θ which are already present in the ‘seed’ model, we introduce two additional parameters $\theta^+$ and θ⁻. Next, we promote φ_θ to an upgraded architecture

$$\begin{align} \widetilde{\phi}_\theta \left( t \right) = \begin{cases} \phi_\theta \left( \theta^+ t\right) &\text{if } t \unicode{x2A7E} 0 \\ \phi_\theta \left( \theta^- t \right) &\text{if } t \lt 0 \end{cases} \, .\end{align} \tag{ 3.14 }$$

The choice of which parameter $\theta^{\pm}$ describes the behavior at t = 0 is immaterial; for concreteness, we have chosen $\theta^+$. By construction, we see that the architecture is independent of θ⁻ at positive times and independent of $\theta^+$ at negative times, so these parameters satisfy the first requirement of our parameter splitting mechanism.

We must also specify the form of the joint density $\widetilde{P} ( \theta^+ , \theta^- , \theta )$ over all of the parameters. Here there is a great deal of freedom. If there are N parameters θ in the seed model, choose any density $p ( \theta^+ , \theta )$ on N + 1 parameters. Then we define

$$\begin{align} \widetilde{P} \left( \theta^+ , \theta^- , \theta \right) = p \left( \theta^+ , \theta \right) \cdot p \left( - \theta^- , \theta \right) \, .\end{align} \tag{ 3.15 }$$

For instance, if one wishes to choose the new density so that its dependence on the seed parameters θ is unchanged, one could take a factorized form

$$\begin{align} p \left( \theta^+ , \theta \right) = \sqrt{P \left( \theta \right) } \cdot f \left( \theta^+ \right) \, ,\end{align} \tag{ 3.16 }$$

for some function f of one variable, where $P ( \theta )$ is the parameter density of the original model. In this case, we see that the density for the new NN-QM takes the simple form $\widetilde{P} ( \theta^+ , \theta^- , \theta ) = P ( \theta ) f ( \theta^+ ) f ( - \theta^- )$. However, this special choice is not required, so we will continue to work with the general density (3.15).

We have engineered the new architecture (3.14) and density (3.15) so that they admit a parameter absorption trick of the form described around equations (3.9) and (3.10). To see this, consider a positive time t > 0 and compare $\widetilde{\phi}_\theta ( t ) = \phi_\theta ( \theta^+ t )$ to its time-reversal $\widetilde{\phi}_\theta ( - t ) = \phi_\theta ( - \theta^- t )$. Clearly the two expressions agree if one identifies $\theta^+ = - \theta^-$. Thus we see that the time-reversal of the input t can be absorbed into a redefinition of the parameters. To make this very explicit, let us define $\theta^{^{\prime} -} = - \theta^-$ and carefully show the parameter absorption mechanism, writing $\widetilde{\phi}_{\theta^+, \theta^-, \theta}$ so that the dependence on all of the parameters is apparent. Then one has

$$\begin{align} \widetilde{\phi}_{\theta^+ , \theta^- , \theta } \left( - t \right) = \phi_\theta \left( - \theta^- t \right) = \phi_\theta \left( \theta^{^{\prime} -} t \right) = \widetilde{\phi}_{\theta^{^{\prime} -} , \theta^- , \theta } \left( t \right) \, .\end{align} \tag{ 3.17 }$$

The point is that, in the final expression of (3.17), the parameter $\theta^+$ which usually sits in the first slot of $\widetilde{\phi}_{\theta^+, \theta^-, \theta}$ has been replaced with $\theta^{^{\prime} -} = - \theta^-$. This procedure allows us to perform the parameter absorption required by condition (3.9) in our splitting mechanism. The second condition (3.10) is ensured by our choice (3.15) for the density. Identifying

$$\begin{align} P_+ \left( \theta^+ , \theta \right) = p \left( \theta^+ , \theta \right) \, , \quad P_- \left( \theta^- , \theta \right) = p \left( - \theta^- , \theta \right) \, ,\end{align} \tag{ 3.18 }$$

to facilitate comparison with the notation of (3.10), clearly one has

$$\begin{align} P_+ \left( - \theta^+ , \theta \right) = P_- \left( \theta^+ , \theta \right) \, , \quad P_- \left( - \theta^- , \theta \right) = P_+ \left( \theta^- , \theta \right) \, ,\end{align} \tag{ 3.19 }$$

which is the other condition required in the parameter splitting mechanism.

By the general argument of the preceding subsection, we therefore conclude that the architecture (3.14) and density (3.15) define a reflection-positive NN-QM given any input network φ_θ. Although this conclusion follows from the results above, it may be instructive to explicitly present the proof of RP for this model to see how the pieces fit together. Given positive times $t_1 , \ldots, t_n \gt 0$ and a function F of n variables, we have

$$\begin{align} \langle F \left( T F \right)^\ast \rangle & = \int \mathrm d \theta \, \left( \int\mathrm d \theta^+ \, p \left( \theta^+ , \theta \right) \, F \left( \widetilde{\phi}_\theta \left( t_1 \right) , \ldots , \widetilde{\phi}_\theta \left( t_n \right) \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm d \theta^- \, p \left( - \theta^- , \theta \right) \, F \left( \widetilde{\phi}_\theta \left( - t_1 \right) , \ldots , \widetilde{\phi}_\theta \left( - t_n \right) \right) \right)^\ast \nonumber \\ & = \int \mathrm d \theta \, \left( \int\mathrm d \theta^+ \, p \left( \theta^+ , \theta \right) \, F \left( \widetilde{\phi}_\theta \left( t_1 \right) , \ldots , \widetilde{\phi}_\theta \left( t_n \right) \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm d \theta^{^{\prime} -} \, p \left( \theta^{^{\prime} -} , \theta \right) \, F \left( \widetilde{\phi}_{\theta^{^{\prime}}} \left( t_1 \right) , \ldots , \widetilde{\phi}_{\theta^{^{\prime}}} \left( t_n \right) \right) \right)^\ast \nonumber \\ & = \int \mathrm d \theta \, \left| \int\mathrm d \theta^+ \, p \left( \theta^+ , \theta \right) \, F \left( \widetilde{\phi}_\theta \left( t_1 \right) , \ldots , \widetilde{\phi}_\theta ( t_n \right) ) \right|^2 \nonumber \\ &\unicode{x2A7E} 0 \, .\end{align} \tag{ 3.20 }$$

The key step is the redefinition of parameters to $\theta^{^{\prime} -} = - \theta^-$ and the change of variables²⁰ in the integral, which upon recognizing that $\theta^+$ and $\theta^{^{\prime} -}$ are dummy variables which are integrated over, allows us to conclude that the two integrals in parentheses are equal and thus the expectation value $\langle F \left( T F \right)^\ast \rangle$ is the integral of a perfect square. This completes the proof of RP for this class of NN-QM models.

We emphasize that, although this parameter splitting mechanism gives one way of engineering RP, it is not ‘generic’ as there are RP processes that do not achieve RP in this way. As we will discuss further below, the procedure for obtaining RP about t = 0 described here does not preserve translation invariance and thus violates (E1). This construction should therefore not be viewed as a complete or universal method for achieving RP, but rather as one new example of a way to guarantee this property, which is interesting since methods of ensuring RP are quite rare.

Translation invariance  and the inevitability of non-analyticity

The construction of $\tilde{\phi}_\theta$ in equation (3.14), with an appropriate parameter density (3.15), produces a reflection positive NN-QM given any input NN φ_θ. However, the resulting model is not time-translation invariant, even if the original network φ_θ enjoys translational symmetry. This is manifest from the form of the architecture (3.14), which explicitly singles out the time t = 0 as special. For example, a correlation function like $\langle \widetilde{\phi}_\theta ( - 1 ) \widetilde{\phi}_\theta ( 1 ) \rangle$ is sensitive to the statistics of both $\theta^+$ and θ⁻, while shifting $t \to t + 2$ gives the correlator $\langle \widetilde{\phi}_\theta ( 1 ) \widetilde{\phi}_\theta ( 3 ) \rangle$, which now only involves $\theta^+$. A priori there is no reason for the shifted and un-shifted correlators to agree, and in general they will not match.

This brings us to an important point about NN-QM systems which satisfy all of the OS axioms. In the preceding discussion we have been focused on RP in the form (3.1), which involves reflections about the point t = 0. However, for a quantum system which enjoys both time-translation invariance and RP, the model must exhibit positivity with respect to reflections about any time t = a. This is a much stronger condition, and not one which can be easily engineered using parameter splitting.

To better understand the tension between the parameter splitting mechanism for RP and translation invariance, let us think about how one might try to engineer RP about two points t₁, t₂ using the above construction. Naively, one would now need to introduce four extra parameters $\theta^{1 \pm}$ and $\theta^{2 \pm}$ with the properties that $\phi_\theta ( t )$ is independent of $\theta^{1 +}$ for $t \lt t_1$, independent of $\theta^{1 -}$ for $t \gt t_1$, independent of $\theta^{2 +}$ for $t \lt t_2$, and independent of $\theta^{2 -}$ for $t \gt t_2$. If one proceeds in this way, one would be naturally led to introduce an uncountable infinity of parameters $\theta^{a \pm}$ for $a \in \mathbb{R}$, which seems excessive.

Another way to think of this tension is to return to our intuitive picture of the parameter splitting mechanism as a ‘decoupling’ between the parameter-space degrees of freedom $\theta^+$ which control the positive-time behavior of the network and the parameters θ⁻ which determine the behavior for negative times. In order to engineer RP about any point t = a via parameter splitting, it would therefore seem that we would need the behavior of the NN at any time t to admit a decoupling from its behavior at other times $t^{^{\prime}} \neq t$. In particular, it appears that the NN outputs at time t should never be determined from the outputs in a neighborhood of t, since this would preclude any such decoupling. This can be made precise in the following result.

Proposition 3.1. Let φ_θ be a NN-QM which is reflection-positive with respect to all times. Assume that the mechanism by which φ_θ achieves RP about each time t is by parameter splitting. Then φ_θ cannot be an analytic function of t at any point in its domain.

Proof. Suppose by way of contradiction that φ_θ is an analytic function of time near a point t = a. This means that, within some finite radius of convergence $| t - a| \lt R$, the outputs of the NN are described by a convergent Taylor series expansion

$$\begin{align} \phi_\theta \left( t \right) = \sum_{n = 0}^{\infty} c_n \left( \theta \right) t^n \, .\end{align} \tag{ 3.21 }$$

Note that φ_θ is still a random function, so different draws of the NN will give different analytic functions near t = a, with different Taylor coefficients that are determined by the values of the parameters. We indicate this dependence by writing the series coefficients as $c_n ( \theta )$. Equation (3.21) defines what we mean by analyticity for NNs.

By assumption, φ_θ achieves RP about t = a by a parameter splitting mechanism. We therefore wish to identify two sets of parameters $\{\theta^+ \}$, $\{\theta^- \}$ such that $\phi_\theta ( t )$ is independent of $\theta^+$ for t < a and $x_\theta ( t )$ is independent of θ⁻ for t > a.

However, for each fixed parameter $\bar{\theta}$, there are only two possibilities: either none of the Taylor coefficients depend on $\bar{\theta}$, or at least one of the c_n depends on this $\bar{\theta}$.

• (a)  
  In the first case, the behavior of φ_θ is independent of this parameter $\bar{\theta}$ for both t < a and t > a, in which case $\bar{\theta}$ cannot be used for parameter splitting. This is because, if φ_θ is independent of $\bar{\theta}$ on both sides of t = a, then (3.9) can never hold.
• (b)  
  In the second case, where at least one of the c_n depends on $\bar{\theta}$, then the behavior of φ_θ depends on $\bar{\theta}$ for both t < a and t > a because the same convergent Taylor series expansion involving the same c_n is used in both regions. We can see this directly by computing a quantity that depends on $\bar{\theta}$ using only $\phi_\theta ( t )$ for t > a, and again using only data for t < a. For instance, one could compute any of the c_n via the formula

  $$\begin{align} c_n \left( \theta \right) = \frac{1}{n!} \frac{\mathrm d^n}{\mathrm d t^n} \phi_\theta \left( t \right) \Big\vert_{t = a} \, ,\end{align} \tag{ 3.22 }$$

  and the derivatives (3.22) could be computed using one-sided limits which only require knowing the behavior of $\phi_\theta ( t )$ for t > a (or, equivalently, only knowing the function for t < a). This shows explicitly that, if c_n depends on some $\bar{\theta}$, then $\phi_\theta ( t )$ depends on $\bar{\theta}$ on both sides of the t = a point. Thus any parameter $\bar{\theta}$ on which at least one of the c_n depends also cannot be used for splitting mechanism.

This argument applies to any parameter, so we conclude that there do not exist any parameters $\theta^+$ and θ⁻ that can be used for the parameter splitting mechanism. This contradicts that φ_θ achieves RP about t = a by parameter splitting. □

Proposition 3.1. may seem disturbing if one is accustomed to thinking of NNs with analytic architectures like the cosine in (2.10), or of conventional NNs such as a feedforward network φ_θ with ReLU activation functions, since this φ_θ is an analytic function of its arguments away from a set of measure zero corresponding to the kink discontinuities of the functions $\mathrm{ReLU} ( z )$ at z = 0. However, in the rigorous formulation of the Euclidean quantum mechanics path integral, one may also recall that the path integral measure is supported on paths which are differentiable nowhere. This construction proceeds by combining the naive path integral measure $\mathcal{D} x$ with the exponentiated kinetic part of the action $\mathrm e^{-S_{\text{kin}}}$, where $S_{\text{kin}} \sim \int \mathrm dt \, \dot{x}^2$, in a particular discretization scheme. Taking the continuum limit of this discretization, one can prove that the combination $\mathrm e^{-S_{\text{kin}}} \, \mathcal{D} x$ converges to the Wiener measure on continuous functions. Functions which are differentiable at even a single point form a set of measure zero with respect to the Wiener measure; the full-measure contribution comes from nowhere-differentiable paths, such as typical sample paths of Brownian motion. From this perspective, it is not so surprising that a reasonable NN-QM system might involve architectures which are analytic nowhere (or differentiable nowhere, which is a stronger condition).

Let us also point out that the inclusion of nowhere-differentiable paths in the standard quantum mechanics path integral is not merely a mathematical curiosity, but has physical implications. For instance, in appendix we review the standard argument that non-trivial commutators such as $[ \hat{x} , \hat{p} ] = i \hbar$ are a direct consequence of the inclusion of non-differentiable trajectories in the path integral. In particular, if one studies a NN-QM system $x ( t ) = \phi_\theta ( t )$ whose architecture is a continuously differentiable function of t at each point, then it would necessarily follow that $[ \hat{x} , \hat{\dot{x}} ] = 0$. Assuming that $\dot{x}$ is related to the canonical momentum, this contradicts the fundamental commutation relations that we expect for a quantum system. Therefore, if we would like our NN-QM to reproduce basic features of quantum mechanics such as commutation relations and the associated uncertainty principles, nowhere-differentiable architectures are essential. However, as we will see in section 4, such nowhere-differentiable paths can be approximated numerically using piecewise-differentiable paths with kink discontinuities at a finite but large collection of points.

We now turn to a Markov process description of the quantum system, where a Markov process y_t could arise as the ‘output-space’ description of a NN with $y_t = x_t$, or in the parameter space perspective if $y_t = \theta_t$. The latter case achieves Markov $x_{\theta_t}$ if the architecture is injective, but we emphasize the former output space viewpoint in this section for simplicity. In this perspective, x_t defines joint probability distributions (2.8) over outputs, rather than the ‘parameter-space’ description in terms of a NN architecture and probability density over parameters θ. This has the advantage of being ‘representation-independent’ in the sense that one need not choose a particular parameterization involving specific θ. As a result, we can present a mechanism for reflection-positivity, somewhat analogous to parameter splitting, which is phrased directly as a sufficient condition on the process x_t and without referring to its parameterization. Furthermore, we will be able to easily implement both RP and the OS axiom (E1) of symmetry and stationarity at once.

To motivation this condition, suppose that we would like to restrict attention to SPs with nowhere-differentiable sample paths, motivated by the importance of such paths for commutation relations which we mentioned at the end of the last subsection. Famous examples of such processes include Brownian motion, Ornstein–Uhlenbeck processes, and any diffusion process obeying a SDE

$$\begin{align} \mathrm d x_t = \mu \left( x_t \right) \, \mathrm dt + \sigma \left( x_t \right) \,\mathrm d W_t \, ,\end{align} \tag{ 3.23 }$$

where W_t is a Wiener path. If $\sigma ( x_t ) \neq 0$, any solution to an equation of the form (3.23) is nowhere-differentiable almost surely. All of these processes share the Markov property, which is the statement that ‘the future is independent of the past, given the present.’ This is clear given an expression (3.23) for the infinitesimal increments $\mathrm d x_t$, which manifestly depend only on the present value of x_t but not on its values at earlier times.

More precisely, a SP x_t is said to be Markov if, given times $t_1 \lt t_2 \lt \ldots \lt t_n \lt t_{n+1}$, one has

$$\begin{align} \mathbb{P} \left( x \left( t_{n+1} \right) \mid x \left( t_1 \right) = x_1 \, , \, x \left( t_2 \right) = x_2 \, , \, \ldots \, , x \left( t_n \right) = x_n \right) = \mathbb{P} \left( x \left( t_{n+1} \right) \mid x \left( t_n \right) = x_n \right) \, ,\end{align} \tag{ 3.24 }$$

where a vertical bar indicates a conditional probability distribution. Said differently, conditioning on a collection of observations $x ( t_i ) = x_i$ of the SP at several past points gives no additional information about the future evolution of the SP than conditioning on the single, latest observation $x ( t_n ) = x_n$.

Intuitively, the Markov property plays a similar role as the parameter splitting mechanism in NN-QM. If a NN-QM $\phi_\theta ( t )$ enjoys parameter splitting with respect to some time t_n, this means that there are parameters θ⁻ which do not affect the outputs of the network for times $t \gt t_n$, although they are relevant for the model’s behavior at earlier times $t \lt t_n$. These parameters θ⁻ are similar to the outputs $x ( t_1 )$, $\ldots$, $x ( t_{n-1} )$ in the Markov condition (3.24). Given the value of $x ( t_n )$, which is roughly analogous to the parameters θ⁰ in the splitting mechanism, future outputs like $x ( t_{n+1} )$ are independent of the behavior of the SP for $t \lt t_n$, just as the future behavior of a NN with parameter splitting is independent of the parameters θ⁻. As a simple example, a solution to a SDE—presented as a NN via its KKL decomposition—is Markov in the ‘output-space’ sense but does not enjoy the ‘parameter-space’ Markov property, since each of the random coefficients θ^k multiplies a function e_k which is generically non-zero at both positive and negative times, and thus the parameters do not split. On the other hand, the extreme case of a network $\phi_\theta ( t ) = \theta_t$, with an uncountable infinity of independent parameters θ_t for each time t, enjoys both the output-space Markov property and the parameter-space Markov property.

Because parameter splitting implies RP, and the Markov property is conceptually similar to parameter splitting, one might ask whether the Markov property is related to RP. The answer is affirmative, and this connection has been known for some time. In fact, the condition of RP was originally introduced in [14] as a weakening of the Markov condition used by Nelson [34]. A proof that the Markov property implies RP, in the general d-dimensional case, can be found in the chapter ‘Probability Theory and Euclidean Field Theory’ by Nelson in the lecture notes [35], and a similar two-line argument (using slightly different language) appears in proposition 1.5 of [36]. Hence it is natural to think of RP as a generalization of the Markov property [37].

Although this result is well-understood, even for d-dimensional Euclidean field theories, we find it useful to give an explicit proof in the d = 1 case of quantum mechanics which we focus on in this work. By virtue of restricting to a one-dimensional spacetime, we are able to give an argument which is more elementary than the one for general d, and which makes the resemblance to the parameter splitting mechanism clear. For a proof of the analogous result for discrete time SPs, see proposition 2.2.5 of [38].

Theorem 3.1 Every symmetric Markov process is reflection positive.

Proof. Let x_t be Markov and symmetric, and fix a collection of positive times $0 \lt t_1 \lt \ldots \lt t_n$. Consider the joint probability distribution

$$\begin{align} \mathbb{P} = \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) , x \left( 0 \right) , x \left( t_1 \right) , \ldots , x \left( t_n \right) \right) \, ,\end{align} \tag{ 3.25 }$$

where we have added the time-zero value x(0) (which is not in the list of t_i and $-t_i$). The quantity $\mathbb{P}$ can be written in conditional form as

$$\begin{align} \mathbb{P} & = \mathbb{P} \left( x \left( t_1 \right) , \ldots , x \left( t_n \right) \mid x \left( - t_n \right) , \ldots, x \left( - t_1 \right) , x \left( 0 \right) \right) \cdot \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) , x \left( 0 \right) \right) \, \nonumber \\ & = \mathbb{P} \left( x \left( t_1 \right) , \ldots , x \left( t_n \right) \mid x \left( 0 \right) \right) \cdot \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) , x \left( 0 \right) \right) \, ,\end{align} \tag{ 3.26 }$$

where in the second line we have used the Markov property. Similarly, one can write

$$\begin{align} \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) , x \left( 0 \right) \right) = \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) \mid x \left( 0 \right) \right) \cdot \mathbb{P} \left( x \left( 0 \right) \right) \, ,\end{align} \tag{ 3.27 }$$

so that

$$\begin{align} \mathbb{P} = \mathbb{P} \left( x \left( t_1 \right) , \ldots , x \left( t_n \right) \mid x \left( 0 \right) \right) \cdot \mathbb{P} \left( x \left( - t_n \right) , \ldots, x \left( - t_1 \right) \mid x \left( 0 \right) \right) \cdot \mathbb{P} \left( x \left( 0 \right) \right) \, .\end{align} \tag{ 3.28 }$$

Now given any appropriate function $F \left( x ( t_1 ) , \ldots , x ( t_n ) \right)$ that depends on the values of x_t at positive times, one has

$$\begin{align} \langle F \left( T F \right)^\ast \rangle & = \int \mathbb{P} \left( x_{-n} , \ldots, x_{-1} , x_0 , x_1 , \ldots , x_n \right) \,\mathrm d x_{-n} \ldots\mathrm d x_{-1} \,\mathrm d x_0 \,\mathrm dx_{1} \ldots\mathrm d x_n \nonumber \\ &\quad \quad \cdot F \left( x_1 , \ldots , x_n \right) \left( F \left( x_{-1} , \ldots , x_{-n} \right) \right)^\ast \, ,\end{align} \tag{ 3.29 }$$

where T is the time-reversal operator and where we use the shorthand notation $x_i = x ( t_i )$, $x_{-i} = x ( - t_{i} )$, $\mathrm d x_{-i} =\mathrm d x ( - t_i )$, $\mathrm d x_i =\mathrm d x ( t_i )$. Using the decomposition (3.28), this is

$$\begin{align} \langle F \left( T F \right)^\ast \rangle & = \int \mathbb{P} \left( x_1 , \ldots , x_n \mid x_0 \right) \cdot \mathbb{P} \left( x_{-n} , \ldots, x_{-1} \mid x_0 \right) \cdot \mathbb{P} \left( x_0 \right) \nonumber \\ &\quad \cdot \,\mathrm d x_{-n} \ldots\mathrm d x_{-1} \,\mathrm d x_0 \,\mathrm dx_{1} \ldots\mathrm d x_n \cdot F \left( x_1 , \ldots , x_n \right) \left( F \left( x_{-1} , \ldots , x_{-n} \right) \right)^\ast \, \nonumber \\ & = \int\mathrm d x_0 \, \mathbb{P} \left( x_0 \right) \, \left( \int\mathrm dx_1 \ldots\mathrm d x_n \, \mathbb{P} \left( x_1 , \ldots , x_n \mid x_0 \right) F \left( x_1 , \ldots , x_n \right) \right) \nonumber \\ &\quad \cdot \left( \int\mathrm dx_{-n} \ldots \mathrm dx_{-1} \, \mathbb{P} \left( x_{-n} , \ldots, x_{-1} \mid x_0 \right) F \left( x_{-1} , \ldots , x_{-n} \right) \right)^\ast \, .\end{align} \tag{ 3.30 }$$

The integral in the parentheses of the first line in the final equation of (3.30) computes a particular conditional expectation value of the function F in the SP x_t. Similarly, the integral in the final line of (3.30) computes the same conditional expectation value in the SP x_−t. Since x_t is symmetric, the processes x_t and x_−t are equivalent and these expectations are equal. We therefore conclude that

$$\begin{align} \langle F \left( T F \right)^\ast \rangle = \int\mathrm d x_0 \, \mathbb{P} \left( x_0 \right) \, \left| \int\mathrm dx_1 \ldots\mathrm d x_n \, \mathbb{P} \left( x_1 , \ldots , x_n \mid x_0 \right) F \left( x_1 , \ldots , x_n \right) \right|^2 \, ,\end{align} \tag{ 3.31 }$$

which is manifestly non-negative. □

Comparing, for instance, the penultimate lines of equations (3.11) and (3.30) highlights the technical similarities between the parameter splitting and Markov mechanisms for achieving RP. In both cases, one factorizes the expectation value $\langle F ( T F )^\ast \rangle$ into an integral of a product of two factors, one of which is another expectation value associated with only positive times and the other of which concerns only negative times. Using an appropriate symmetry—arising from either the parameter absorption trick, in the splitting mechanism, or from the explicit assumption that x_t and x_−t are equivalent, in the Markov argument—one then observes that these two factors are equal, which means that $\langle F ( T F )^\ast \rangle$ is the integral of a perfect square, and is hence positive.

However, a notable difference between the two mechanisms is the ease or difficulty of combining RP with translation invariance. We saw in section 3.1 that incorporating time translation symmetry into the parameter splitting mechanism seems tricky, and naively appears to require introducing an uncountable infinity of parameters. However, from the ‘output-space’ perspective, combining these two features is easy. Indeed, we immediately see that every symmetric, stationary Markov process is positive with respect to reflections about any time t = a, and thus satisfies both of the OS axioms (E1) and (E4). Morally, this is because the definition (3.24) of a Markov process makes no reference to any specific time t = 0. Rather, the ‘independence of the future from the past, given the present’ occurs universally for any value of the present time.

First let us discuss the Markov property, which was used to achieve RP in theorem 3.1. Unfortunately, the space of Markov processes is rather ‘fragile’ in the sense that it is not closed under many operations. For instance, a linear combination of Markov processes need not be Markov. Intuitively, this is because knowing the sum of two numbers provides strictly less information than knowing each of the two summands individually. If $z_t = x_t + y_t$ where x_t, y_t are Markov, then knowing the present values of x_t and y_t gives enough information to predict their future behavior, and hence the behavior of their sum z_t. But if one only knows the present value of z_t, it is impossible to recover the precise current values of x_t and y_t since many pairs of numbers have the same sum, and thus the present value of z_t is not sufficient to predict its future behavior.

Similarly, the space of Markov processes is not closed under the operation of applying deterministic functions. As an illustrative example, if W_t is a Wiener process or Brownian path, [39] classified all functions $f : \mathbb{R} \to \mathbb{R}$ for which $x_t = f ( W_t )$ is Markov, and found only four families of possibilities. Again, this is because a general function is not injective, so knowing the value of $f ( W_t )$ is typically insufficient to recover the value of W_t.

However, it was pointed out in [37] that the space of RP processes behaves more nicely with respect to these operations. In particular, one has the following results.

Proposition 4.1. If x_t and y_t are two reflection positive processes, then any linear combination $z_t = a x_t + b y_t$, for $a, b \in \mathbb{R}$, is also reflection positive.

Proposition 4.2. If x_t is reflection positive and $f : \mathbb{R} \to \mathbb{R}$ is a bounded measurable deterministic function, then the process $z_t = f ( x_t )$ is reflection positive.

It is also straightforward to see that the operations of taking linear combinations and/or deterministic functions preserve the properties of symmetry and . These observations suggests a strategy for generating a large class of Euclidean quantum models that admit real-time continuations: begin with some collection of simple examples, such as symmetric stationary Markov processes, and apply a sequence of such operations to produce new systems which still obey the OS axioms (E1) and (E4).

One can generalize proposition 4.2 further by considering random functions of RP processes, such as those defined by a NN, which still preserve RP. This is the main result of this section, and is formalized in the following theorem.

Theorem 4.1. Let x_t be a reflection-positive SP and consider a family of measurable random functions realized by a NN architecture φ_θ with parameter density $P ( \theta )$. Then the SP

$$\begin{align} y_t = \phi_\theta \left( x_t \right)\end{align} \tag{ 4.1 }$$

is also reflection-positive.

Note that y_t has two sources of randomness: to generate a sample path of y_t, we first generate an entire realization of the process x_t, and we perform a random draw of the parameters θ to obtain an instance of the function φ_θ, then apply φ_θ to x_t pointwise.

Proof. Fix a collection of positive times $0 \lt t_1 \lt \ldots \lt t_n$. Following the notation in the proof of theorem 4.1, we write $y_i = y ( t_i )$, $y_{-i} = y ( - t_i )$, and so on. Let $F ( y_1 , \ldots, y_n )$ be a bounded measurable function which depends only on the values of $y_i = y ( t_i )$ at these positive times. Then consider the expectation value

$$\begin{align} \langle F \left( T F \right)^\ast \rangle_y & = \int \mathrm d \theta \, P \left( \theta \right) \int\mathrm d x_{-n} \ldots\mathrm d x_{-1} \,\mathrm dx_1 \, \ldots \,\mathrm d x_n \, \mathbb{P} \left( x_{-n} , \ldots, x_n \right) \nonumber \\ &\quad F \left( \phi_\theta \left( x_1 \right) , \ldots, \phi_\theta \left( x_n \right) \right) \left( F \left( \phi_\theta \left( x_{-1} \right) , \ldots, \phi_\theta \left( x_{-n} \right) \right) \right)^\ast \, .\end{align} \tag{ 4.2 }$$

We use the notation $\langle \, \cdot \, \rangle_y$ to emphasize that this is an expectation value taken with respect to the process $y_t = \phi_\theta ( x_t )$. Now define the function $\widetilde{F} = F \circ \phi_\theta$ of n variables by

$$\begin{align} \widetilde{F} \left( x_1 , \ldots , x_n \right) = F \left( \phi_\theta \left( x_1 \right) , \ldots, \phi_\theta \left( x_n \right) \right) \, ,\end{align} \tag{ 4.3 }$$

so that

$$\begin{align} \langle F \left( T F \right)^\ast \rangle_y & = \int \mathrm d \theta \, P \left( \theta \right) \int\mathrm d x_{-n} \ldots\mathrm d x_{-1} \,\mathrm dx_1 \, \ldots \,\mathrm d x_n \, \mathbb{P} \left( x_{-n} , \ldots, x_n \right) \nonumber \\ &\quad\widetilde{F} \left( x_1 , \ldots , x_n \right) \left( \widetilde{F} \left( x_{-1} , \ldots , x_{-n} \right) \right)^\ast \nonumber \\ & = \int \mathrm d \theta P \left( \theta \right) \, \langle \widetilde{F} \left( T \widetilde{F} \right)^\ast \rangle_x \, .\end{align} \tag{ 4.4 }$$

We have introduced the symbol $\langle \, \cdot \, \rangle_x$ to refer to the expectation value taken with respect to the SP x_t. For each fixed choice of the parameters θ, the function $\widetilde{F}$ is a bounded measurable function of n variables, and thus satisfies the technical condition on the function appearing in the statement (3.1) of RP. Since we have assumed the SP x_t is reflection positive, one has

$$\begin{align} \langle \widetilde{F} \left( T \widetilde{F} \right)^\ast \rangle_x \unicode{x2A7E} 0 \, .\end{align} \tag{ 4.5 }$$

We therefore see that $\langle F ( T F )^\ast \rangle_y$ is the integral of a non-negative definite quantity, multiplied by a non-negative parameter density $P ( \theta )$, so we conclude that

$$\begin{align} \langle F \left( T F \right)^\ast \rangle_y \unicode{x2A7E} 0 \, ,\end{align} \tag{ 4.6 }$$

and hence y_t is reflection positive. □

Theorem 4.1 motivates the following construction. Let $y_t^{(1)}$, $\ldots$, $y_t^{(n)}$ be a finite set of reflection-positive processes (for instance, symmetric Markov processes). If the $y_t^{(i)}$ are MQMs, we may realize each of them as NNs by presenting them using the KKL decomposition. We then view these processes as inputs to a new NN. Choose a width D₁ for the next layer and define a collection of D₁ processes $z^{(1, i)}_t$ by

$$\begin{align} z^{\left(1, i\right)}_t = \phi_\theta^{\left(1, i\right)} \left( \sum_{j = 1}^{n} c_{j}^{\left(1, i\right)} y_t^{\left(j\right)} \right) \, , \quad i = 1 , \ldots, D_1 \, ,\end{align} \tag{ 4.7 }$$

where $\phi_\theta^{(1, i)}$ are NNs (or general random functions) and $c_{j}^{(1, i)}$ are real coefficients, which can be chosen either to be constants or random variables with some probability distribution. We then continue: choosing another width D₂, we define

$$\begin{align} z^{\left(2, i\right)} = \phi_\theta^{\left(2, i\right)} \left( \sum_{j = 1}^{D_1} c_{j}^{\left(2, i\right)} z_t^{\left(j\right)} \right) \, , \quad i = 1 , \ldots, D_2 \, ,\end{align} \tag{ 4.8 }$$

and so on, up to a collection of processes $z^{(N, i)}$, $i = 1 , \ldots, D_N$, for a positive integer N. Finally, we construct the output layer

$$\begin{align} x_t = \phi_\theta^{\left(x\right)} \left( \sum_{j = 1}^{D_N} c_{j}^{\left(N\right)} z_t^{\left(N, j\right)} \right) \, .\end{align} \tag{ 4.9 }$$

We refer to any quantum system x_t obtained by a sequence of steps of this form as a deep NN quantum mechanics (deep NN-QM). For a fixed choice of architecture of the type described above, and fixed input processes, changing the parameter density $P ( \theta )$ substantially modifies the properties of the output process; the choice of density is part of the definition of the deep NN-QM theory. Our theorem guarantees that RP is preserved at each layer of this network, so the resulting SP x_t is reflection positive. If all of the input processes $y_t^{(i)}$ are symmetric and stationary, so that they satisfy the OS axiom (E1), then x_t enjoys this property as well. This is because the NN architecture and parameters are independent of time, and the function defining x_t is obtained by applying operations to the inputs $y_t^{(i)}$ pointwise (at each value of time), so this construction does not break any symmetries involving transformation of the argument t, such as $t \to t + a$ or $t \to - t$²¹.

Although it is not our primary focus here, let us remark that certain choices of deep NN-QM architectures also preserve the mixing property, in the sense that if each of the $y_t^{(i)}$ obeys the cluster decomposition axiom (E3), then x_t also inherits this property. This is not, however, true in general: a counter-example is a network which simply shifts the input process by a random offset, $y_t = x_t + \theta$ with θ drawn randomly and separately for each sample path x_t. As is known in ergodic theory, such an ensemble of mixing processes with random shifts may fail to be mixing (this is in contrast with the application of a deterministic measurable function, which always preserves the mixing property).

Therefore, given any collection of input process $y_t^{(i)}$ which satisfy all of the OS axioms, the deep NN-QM construction generically produces a new quantum system that also satisfies (E0), (E1), (E2), and (E4), and for certain choices also (E3); in the latter case the deep NN-QM hence gives rise to a well-defined real-time quantum mechanical theory.

The quantum systems produced by this deep NN-QM procedure are non-trivial, even if the input processes $y_t^{(i)}$ are quite simple. For instance, suppose that all of the $y_t^{(i)}$ are Ornstein–Uhlenbeck processes, which are associated with quantum harmonic oscillators with some mass parameter. Then in addition to being mixing, symmetric, stationary, and Markov—and consequently satisfying all of the OS axioms—the systems $y_t^{(i)}$ are Gaussian, and are in some sense trivial since they are analogous to free quantum field theories. However, the output x_t of a generic deep NN-QM built from these Markov, Gaussian input processes will be neither Markov nor Gaussian, although such an x_t is symmetric, stationary, and RP. The non-Gaussianity of x_t is clear since a non-linear transformation of a Gaussian process is no longer Gaussian, and the NNs $\phi_\theta^{(i, j)}$ used in this construction are nonlinear in general. In particular, the higher-order cumulants of a process obtained from applying a generic deep NN-QM to Gaussian inputs will be non-vanishing, which signals that they are interacting quantum mechanical theories.

In this section, we will present some numerical results for quantum observables obtained by simulating SPs on a computer. The simplest model which can be studied in this way is the quantum harmonic oscillator, which is associated with the statistics of the Ornstein–Uhlenbeck process and is discussed in section 4.2.1. We then use OU processes as inputs for a deep NN-QM construction in section 4.2.2. As we have mentioned, since the Ornstein–Uhlenbeck process is symmetric, stationary, and Markov, any such deep NN-QM constructed from OU building blocks is symmetric, stationary, and reflection positive. We will choose examples of deep NN-QM architectures for which the output process is also mixing, which gives rise to families of Euclidean quantum models obeying all of the OS axioms and which can thus be continued to real time.

4.2.1. Ornstein–Uhlenbeck process

The Ornstein–Uhlenbeck process [40] is a famous and well-studied SP, which enjoys many desirable features that make it amenable to theoretical analysis, while still exhibiting sufficient richness as to be interesting. One can think of the OU process as a mean-reverting random walk, in the sense that the trajectory x_t of such a process exhibits Brownian fluctuations in addition to a ‘restoring force’ which acts to return the position to its equilibrium value $x_t = 0$. This two effects are modeled mathematically via the two terms on the right side of the SDE

$$\begin{align} \mathrm dx_t = - \theta x_t \, \mathrm dt + \sigma \,\mathrm d W_t \, ,\end{align} \tag{ 4.10 }$$

where W_t is a Wiener path. Here θ measures the strength of the deterministic ‘drift’ part of the SDE, which leads to the mean-reverting behavior, and the ‘noise’ parameter σ describes the size of the random fluctuations. This is an example of the general form (3.23) of the SDE which defines an Itô diffusion. As we have remarked, the solution of any such SDE is a Markov process which, assuming symmetry, is RP.

An OU process with a fixed initial condition, say $x ( 0 ) = 0$, is not stationary. However, the Ornstein–Uhlenbeck process admits a stationary distribution, in the sense that if the initial condition $x_0 = x ( 0 )$ is drawn from a particular probability distribution $\mathbb{P}_{s} ( x_0 )$, then the resulting SP is stationary. This distribution is given by

$$\begin{align} \mathbb{P}_{s} \left( x_0 \right) = \sqrt{\frac{\theta}{\pi \sigma^2} } \exp \left( - \frac{\theta x_0^2}{\sigma^2} \right) \, ,\end{align} \tag{ 4.11 }$$

which is a Gaussian distribution with mean zero and variance $\frac{\sigma^2}{2 \theta}$. When we speak of ‘the’ Ornstein–Uhlenbeck process, we always mean the version of this process where the initial position of the process is drawn from the stationary distribution (4.11). In this case, the OU process is stationary, symmetric, and Markov, and it also turns out to be Gaussian²².

The statistics of the Ornstein–Uhlenbeck process precisely correspond to the Euclidean-time quantum harmonic oscillator, with a Hamiltonian that has been shifted by a constant so that the ground state energy is E₀ = 0 rather than $E_0 = \frac{1}{2} \hbar \omega$. As we will review in more detail around equation (4.38), this can be seen by performing a similarity transformation on the Fokker–Planck operator of the Ornstein–Uhlenbeck process [42], bringing it to the form of an effective Hamiltonian

$$\begin{align} H_{\text{OU}} = - \frac{\sigma^2}{2} \frac{\mathrm d^2}{\mathrm dx^2} + \frac{\theta^2}{2 \sigma^2} x^2 - \frac{\theta}{2} \, .\end{align} \tag{ 4.12 }$$

Comparing this to the Euclidean-time Hamiltonian of the quantum harmonic oscillator,

$$\begin{align} H_{\text{SHO}} = - \frac{\hbar^2}{2m} \frac{\mathrm d^2}{\mathrm dx^2} + \frac{1}{2} m \omega^2 x^2 \, ,\end{align} \tag{ 4.13 }$$

we see that the dictionary between the parameters is

$$\begin{align} \theta = \omega \, , \quad \sigma^2 = \frac{\hbar}{m} \, ,\end{align} \tag{ 4.14 }$$

and that the two Hamiltonians are related by

$$\begin{align} H_{\text{OU}} = \frac{1}{\hbar} \left( H_{\text{SHO}} - \frac{1}{2} \hbar \omega \right) \, .\end{align} \tag{ 4.15 }$$

The overall factor of $\frac{1}{\hbar}$ is due to the fact that the quantum mechanical time evolution operator is $\mathrm e^{- H t / \hbar}$, whereas the generator of time translations in the SP has no factor of $\hbar$. We therefore see that the two Hamiltonians H_OU and H_SHO agree, up to this overall scaling and the shift by $\frac{\theta}{2} = \frac{\hbar \omega}{2}$ which we mentioned.

In our numerical analysis of the OU process, we will generate sample paths by solving the SDE (4.10) using the Euler–Maruyama method. This means that we partition a finite interval $[0, T]$ into a collection of points $0 \lt t_1 \lt t_2 \lt \ldots \lt t_N = T$, differing by a time step $t_{i+1} - t_i = \Delta t$, and recursively define the value of $x_i = x_{t_i}$ at these points t_i using

$$\begin{align} x_{i+1} = x_{i} - \theta x_{i} \, \Delta t + \sigma \sqrt{\Delta t } \zeta_i \, ,\end{align} \tag{ 4.16 }$$

where $\zeta_i \sim \mathcal{N} ( 0, 1 )$ are unit normal random variables. The initial condition x₀ is drawn from the stationary distribution (4.11) of the OU process. There are several sources of numerical error in computing observables using a simulation scheme of this form. One is Monte Carlo error: we generate a finite number N of sample paths and compute ensemble averages over quantities in each path, which leads to sampling error for quantities like the two-point function (to be discussed shortly) that scales as $\frac{1}{\sqrt{N}}$. Another source of error is due to the discretization Δt, as the true properties of the continuous OU process are reproduced only as $\Delta t \to 0$. When we compute the spectrum, we will create a histogram of pairs $(x_0 , x_{\Delta t})$ to be explained below, which introduces binning error that also depends on Δt. We will not undertake a systematic analysis of these various sources of error, but in the plots of our numerical results we will include error bars that represent the standard deviation across different epochs, and these error bars give a sense of the magnitude of the combined error from these various sources.

Two-point function

One of the simplest quantities that one can compute numerically for the OU process is the two-point function $G^{(2)} ( t, s )$. The well-known closed-form expression for this correlator is

$$\begin{align} G^{\left(2\right)} \left( t, s \right) = \langle x \left( t \right) x \left( s \right) \rangle = \frac{\sigma^2}{2 \theta} \exp \left( - \theta | t - s | \right) \, .\end{align} \tag{ 4.17 }$$

Because the Ornstein–Uhlenbeck process is Gaussian, all of the higher correlation functions $G^{(n)}$, $n \unicode{x2A7E} 3$, are completely fixed by the two-point function. The fact that the two-point function decomposes into a product of one-point functions at large separation,

$$\begin{align} \lim_{| t - s | \to \infty} \langle x \left( t \right) x \left( s \right) \rangle = \langle x \left( t \right) \rangle \langle x \left( s \right) \rangle = 0 \, ,\end{align} \tag{ 4.18 }$$

is an implication of the mixing property of the OU process, which is physically interpreted as the statement that the quantum harmonic oscillator enjoys cluster decomposition.

As a first test of our numerics, we compare the simulated two-point function for the Ornstein–Uhlenbeck process to the theoretical expression (4.17) in figure 2, finding agreement between the numerical estimate and the theoretical prediction.

Zoom In Zoom Out Reset image size

In view of the spectral representation (2.5) for the two-point function and the fact that the matrix element $\langle 0 \mid x \mid n \rangle$ vanishes for the harmonic oscillator unless n = 1, we can immediately read off the gap between the ground state energy $E_0^{\text{OU}}$ and the first excited state energy $E_1^{\text{OU}}$ from the rate of the exponential decay of this two-point function. In this case, we find $E_1^{\text{OU}} - E_0^{\text{OU}} = \theta$. Using the dictionary (4.14) and the relation $E_n^{\text{OU}} = \frac{1}{\hbar} E_n^{\text{SHO}}$, which follows from (4.15), this gives $E_1^{\text{SHO}} - E_0^{\text{SHO}} = \hbar \omega$, which agrees with the gap to the first excited state in the harmonic oscillator. We will review a different method for extracting spectral data about the quantum oscillator from the OU process below.

Commutator

We have mentioned, in section 3.1 and appendix, that the existence of nowhere-differentiable paths is crucial for non-trivial commutation relations. We would now like to see how this works for the Ornstein–Uhlenbeck realization of the harmonic oscillator.

Since the Ornstein–Uhlenbeck process satisfies a SDE with non-zero noise term, in the continuum limit the sample paths of this process will be differentiable nowhere almost surely. This means that, in a numerical simulation of OU paths using the Euler–Maruyama method with step size Δt, the estimated derivative

$$\begin{align} \frac{\Delta x}{\Delta t} = \frac{x_{t + \Delta t} - x_t}{\Delta t}\end{align} \tag{ 4.19 }$$

diverges as the discretization becomes finer and finer ($\Delta t \to 0$). However, the expression

$$\begin{align} C \left( t \right) = x_t \frac{x_t - x_{t-\Delta t}}{\Delta t} - x_t \frac{x_{t + \Delta t} - x_t}{\Delta t} \, ,\end{align} \tag{ 4.20 }$$

which estimates the commutator $[ \hat{x} ( t ) , \hat{p} ( t ) ]$, should remain finite as the step size is taken small (as we have argued in appendix for the case of the free quantum particle). Since $[\hat{x} ( t ) , \hat{p} ( t ) ] = \hbar$ for quantum mechanics in Euclidean signature, we expect C(t) to be constant in time and proportional to $\hbar$, which scales as σ² according to the dictionary between OU and SHO parameters given in equation (4.14). In our numerical simulations we find that this expectation is indeed borne out; the results are presented in figure 3. The computation of C(t) gives one way of quantitatively measuring the ‘roughness’ of the sample paths of the Ornstein–Uhlenbeck process²³.

Zoom In Zoom Out Reset image size

Uncertainty  relations

One of the most basic features of quantum mechanics is the Heisenberg uncertainty principle (HUP). In terms of the uncertainty ΔA associated with an observable $\hat{A}$,

$$\begin{align} \left( \Delta A \right)^2 = \langle \hat{A}^2 \rangle - \langle \hat{A} \rangle^2 \, ,\end{align} \tag{ 4.21 }$$

the HUP bounds the product of ΔA and the uncertainty ΔB in a second observable $\hat{B}$:

$$\begin{align} \Delta A \, \Delta B \unicode{x2A7E} \frac{1}{2} \left| \left\langle \left[ \hat{A} , \hat{B} \right] \right\rangle \right| \, .\end{align} \tag{ 4.22 }$$

As we have reviewed, in Euclidean signature the position and momentum operators obey $[ \hat{x} , \hat{p} ] = \hbar$. Therefore, one obtains the general lower bound

$$\begin{align} \Delta x \, \Delta p \unicode{x2A7E} \frac{\hbar}{2} \, .\end{align} \tag{ 4.23 }$$

In the ground state of the harmonic oscillator, which corresponds to the stationary distribution of the OU process, this inequality is saturated and one finds $\Delta x \, \Delta p = \frac{\hbar}{2}$.

We would like to verify the saturation of this uncertainty relation numerically using simulations of the Ornstein–Uhlenbeck process. However, to do this, we must first identify a suitable definition of expectation values of powers of the momentum operator $\hat{p}$. When m = 1, we expect that $p = \dot{x}$, and we have already pointed out that the continuum Ornstein–Uhlenbeck process has sample paths which are differentiable nowhere. We therefore cannot use the naive definition of the derivative $\dot{x}$ for computing momentum expectation values, as these derivatives will diverge as the step size Δt is taken to zero.

There are two ways to sidestep this issue, which turn out to be equivalent for the harmonic oscillator. One way is to represent the momentum operator in terms of the derivative $\frac{\partial}{\partial x}$ acting on wavefunctions, as we do in ordinary quantum mechanics, and compute expectation values of this derivative using the stationary distribution of the Ornstein–Uhlenbeck process. The second way, which we follow here, is motivated by the observation that the OU two-point function (4.17) is non-differentiable at the point t = s, but admits a one-sided derivative. That is, if we restrict to t > s, one has

$$\begin{align} \langle x \left( t \right) x \left( s \right) \rangle = \frac{\sigma^2}{2 \theta} \exp \left( - \theta \left( t - s \right) \right) \, .\end{align} \tag{ 4.24 }$$

One may then take derivatives of this expectation value to define momentum correlators. Following this strategy, we define the second moment of the momentum operator as

$$\begin{align} \langle p \left( t \right)^2 \rangle & = \lim_{\substack{h \to 0^+ \\ k \to 0^+}} \left( \frac{\mathrm d}{\mathrm d h} \frac{\mathrm d}{\mathrm d k} \langle x \left( t + h \right) x \left( t - k \right) \rangle \right) \, .\end{align} \tag{ 4.25 }$$

The notation $h \to 0^+$, $k \to 0^+$ mean that the variables $h, k$ are taken to zero from above (i.e. through strictly positive values). This definition ensures that the ordering

$$\begin{align} t - k \lt t \lt t + h\end{align} \tag{ 4.26 }$$

is maintained throughout the limiting procedure, so that we avoid the non-differentiability of the two-point function at coincident points. In numerical simulations, we use a finite-difference approximation to the quantity (4.25):

$$\begin{align} \langle p \left( t \right)^2 \rangle = \left\langle \frac{x \left( t + \Delta t \right) x \left( t - \Delta t \right) - x \left( t + \Delta t \right) x \left( t \right) - x \left( t \right) x \left( t - \Delta t \right) + x\left(t\right)^2}{\left(\Delta t \right)^2} \right\rangle \, .\end{align} \tag{ 4.27 }$$

Similarly, we define

$$\begin{align} \langle p \left( t \right) \rangle = \frac{\mathrm d}{\mathrm dt} \langle x \left( t \right) \rangle \, ,\end{align} \tag{ 4.28 }$$

and as $\langle x ( t ) \rangle = 0$ this quantity vanishes for the true continuum OU process.

The results of the numerical computation of $\left( \Delta x \right)^2$ and $\left( \Delta p \right)^2$, for various values of θ, σ, and the discretization time Δt, are presented in table 1, which confirms that the uncertainty relation is saturated for the harmonic oscillator ground state within the numerical error induced by our simulation scheme.

Table 1. We compare the simulated (‘sim’) values of $\left( \Delta x \right)^2$, $\left( \Delta p \right)^2$ (defined using (4.25) and (4.28)), and their product, to the predicted theoretical (‘th’) values $\left( \Delta x \right)^2 = \frac{\sigma^2}{2 \theta}$, $\left( \Delta p \right)^2 = \frac{\sigma^2 \theta}{2}$, and $(\Delta x)^2\cdot(\Delta p)^2 = \frac{\sigma^4}{4}$. Simulations are performed using 200 epochs with $10\,000$ OU sample path per epoch, where each path consists of 1000 time steps. We average across time steps and epochs, and the errors indicated by ± are the standard deviation across epochs. We contrast two different sizes of the time step used in the Euler–Maruyama method, $\mathrm dt = 0.1$ and $\mathrm dt = 0.01$, and we see that the simulation results become closer to the theoretical values for the smaller choice of $\mathrm dt$ due to reduced discretization error. The results confirm the relation $\left( \Delta x \right)^2 \left( \Delta p \right)^2 = \frac{\hbar^2}{4}$ for the ground state of the harmonic oscillator, where the effective $\hbar$ of the process is determined by σ² as in (4.14).

θ	σ	$\mathrm dt$	$(\Delta x)^2$ (sim)	$(\Delta x)^2$ (th)	$(\Delta p)^2$ (sim)	$(\Delta p)^2$ (th)	$(\Delta x)^2\cdot(\Delta p)^2$ (sim)	$(\Delta x)^2\cdot(\Delta p)^2$ (th)
$\frac{1}{2}$	$\frac{1}{2}$	0.100	$0.2564 \pm 0.0005$	0.2500	$0.0641 \pm 0.0008$	0.0625	$0.0164 \pm 0.0002$	0.0156
$\frac{1}{2}$	$\frac{1}{2}$	0.010	$0.2509 \pm 0.0015$	0.2500	$0.0633 \pm 0.0077$	0.0625	$0.0159 \pm 0.0019$	0.0156
$\frac{1}{2}$	1	0.100	$1.0256 \pm 0.0018$	1.0000	$0.2560 \pm 0.0034$	0.2500	$0.2626 \pm 0.0034$	0.2500
$\frac{1}{2}$	1	0.010	$1.0023 \pm 0.0064$	1.0000	$0.2503 \pm 0.0316$	0.2500	$0.2509 \pm 0.0316$	0.2500
$\frac{1}{2}$	2	0.100	$4.1004 \pm 0.0077$	4.0000	$1.0270 \pm 0.0125$	1.0000	$4.2109 \pm 0.0513$	4.0000
$\frac{1}{2}$	2	0.010	$4.0080 \pm 0.0261$	4.0000	$1.0099 \pm 0.1334$	1.0000	$4.0475 \pm 0.5330$	4.0000
1	$\frac{1}{2}$	0.100	$0.1316 \pm 0.0002$	0.1250	$0.1315 \pm 0.0008$	0.1250	$0.0173 \pm 0.0001$	0.0156
1	$\frac{1}{2}$	0.010	$0.1257 \pm 0.0005$	0.1250	$0.1252 \pm 0.0075$	0.1250	$0.0157 \pm 0.0009$	0.0156
1	1	0.100	$0.5262 \pm 0.0007$	0.5000	$0.5263 \pm 0.0033$	0.5000	$0.2770 \pm 0.0017$	0.2500
1	1	0.010	$0.5023 \pm 0.0023$	0.5000	$0.5073 \pm 0.0307$	0.5000	$0.2548 \pm 0.0155$	0.2500
1	2	0.100	$2.1044 \pm 0.0029$	2.0000	$2.1072 \pm 0.0137$	2.0000	$4.4343 \pm 0.0277$	4.0000
1	2	0.010	$2.0097 \pm 0.0090$	2.0000	$2.0036 \pm 0.1307$	2.0000	$4.0269 \pm 0.2646$	4.0000
2	$\frac{1}{2}$	0.100	$0.0694 \pm 0.0001$	0.0625	$0.2778 \pm 0.0008$	0.2500	$0.0193 \pm 0.0001$	0.0156
2	$\frac{1}{2}$	0.010	$0.0631 \pm 0.0002$	0.0625	$0.2525 \pm 0.0080$	0.2500	$0.0159 \pm 0.0005$	0.0156
2	1	0.100	$0.2777 \pm 0.0003$	0.2500	$1.1114 \pm 0.0032$	1.0000	$0.3087 \pm 0.0009$	0.2500
2	1	0.010	$0.2524 \pm 0.0009$	0.2500	$1.0154 \pm 0.0317$	1.0000	$0.2563 \pm 0.0079$	0.2500
2	2	0.100	$1.1109 \pm 0.0011$	1.0000	$4.4458 \pm 0.0143$	4.0000	$4.9390 \pm 0.0155$	4.0000
2	2	0.010	$1.0098 \pm 0.0032$	1.0000	$4.0375 \pm 0.1261$	4.0000	$4.0772 \pm 0.1263$	4.0000

Spectrum and eigenfunctions

Finally, we turn to the numerical computation of the energy eigenvalues and eigenfunctions of the harmonic oscillator. Because the spectrum and energy eigenstates are properties of the Hamiltonian operator, which is the generator of infinitesimal time translations, one might expect that this data can be extracted from the update rule (4.10) that describes how OU sample paths evolve over a small time step. Indeed this is the case—by a standard argument, which we will briefly review, the generator of a SP which admits a stationary distribution can be transformed into a Schrodinger-type operator via a certain similarity transformation, which can then be diagonalized to obtain the energies. See, for instance, section 4.9 of [42] for a more thorough discussion of this point.

We first mention that, rather than describing the time evolution of a specific sample path x_t via the SDE (4.10), we can equivalently specify a differential equation that models the dynamics of the probability density $\mathbb{P} ( x, t )$ of an ensemble of OU paths at a time t. This equation takes the form

$$\begin{align} \frac{\partial \mathbb{P}}{\partial t} & = \theta \frac{\partial}{\partial x} \left( x \mathbb{P} \right) + \frac{\sigma^2}{2} \frac{\partial^2 \mathbb{P}}{\partial x^2} \nonumber \\ & = \mathcal{L} P \, ,\end{align} \tag{ 4.29 }$$

which is called the Fokker–Planck equation, and

$$\begin{align} \mathcal{L} f = \frac{\sigma^2}{2} \frac{\partial^2 f}{\partial x^2} + \theta \frac{\partial}{\partial x} \left( x f \right)\end{align} \tag{ 4.30 }$$

is called the generator of the process, or the Fokker–Planck operator. This operator $\mathcal{L}$ generates infinitesimal time translations of the probability density $\mathbb{P} ( x, t )$, but this is not the same as the operator that time-evolves the wavefunction $\psi ( x, t )$ in conventional quantum mechanics, since $\mathbb{P} ( x, t ) = \left| \psi ( x, t ) \right|^2$. One way to see this is to note that the Fokker–Planck operator is not self-adjoint with respect to the inner product on functions:

$$\begin{align} \langle f , \mathcal{L} g \rangle & = \int\mathrm dx \, f \left( x \right) \left( \frac{\sigma^2}{2} \frac{\partial^2 g}{\partial x^2} + \theta g + \theta x \frac{\partial g}{\partial x} \right) \nonumber \\ &\neq \int\mathrm dx \, \left( \frac{\sigma^2}{2} \frac{\partial^2 f}{\partial x^2} + \theta f + \theta x \frac{\partial f}{\partial x} \right) g\left(x\right) \nonumber \\ & = \langle \mathcal{L} f , g \rangle \, .\end{align} \tag{ 4.31 }$$

On the other hand, we expect that the Hamiltonian operator should be self-adjoint when acting on wavefunctions. The failure of self-adjointness is due to the final term $\theta x \frac{\partial g}{\partial x}$ in the first line of (4.31); performing an integration by parts on this term does not yield the expression on the second line. The standard technique for removing such a first derivative term is through the introduction of an integrating factor. Given a differential operator

$$\begin{align} D F = a F^{^{\prime\prime}}\left(x\right) + p \left( x \right) F^{^{\prime}} \left( x \right) + q \left( x \right) F \left( x \right) \, ,\end{align} \tag{ 4.32 }$$

for a constant a and functions p(x) and q(x), one can change variables to a new function

$$\begin{align} G \left( x \right) = \exp \left( \frac{1}{2 a} \int_0^{x} p \left( y \right) \,\mathrm dy \right) F \left( x \right) \, ,\end{align} \tag{ 4.33 }$$

in terms of which the differential operator takes the form

$$\begin{align} D F = \exp \left( - \frac{1}{2 a} \int_0^{x} p \left( y \right) \,\mathrm dy \right)\left( a G^{^{\prime\prime}}\left(x\right) + \left( q \left( x \right) - \frac{1}{2} p^{^{\prime}}\left(x\right) - \frac{1}{4a } p \left( x \right)^2 \right) G \left( x \right) \right) \, ,\end{align} \tag{ 4.34 }$$

which has no term proportional to G′.

In the case of the operator $D = \mathcal{L}$ defined in equation (4.30), taking $a = \frac{\sigma^2}{2}$ and $p(x) = \theta x$, this suggests the change of variables

$$\begin{align} \psi \left( x, t \right) = \exp \left( \frac{\theta}{2 \sigma^2} x^2 \right) \mathbb{P} \left( x, t \right) \, ,\end{align} \tag{ 4.35 }$$

so that the Fokker–Planck operator $\mathcal{L}$ acts as

$$\begin{align} \mathcal{L} \mathbb{P} = \exp \left( - \frac{\theta}{2 \sigma^2} x^2 \right) \left( \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} - \frac{\theta^2}{2 \sigma^2} x^2 + \frac{\theta}{2} \right) \psi \left( x, t\right) \, ,\end{align} \tag{ 4.36 }$$

and in terms of which the differential equation (4.29) can be written as

$$\begin{align} \frac{\partial \psi}{\partial t} = \left( \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} - \frac{\theta^2}{2 \sigma^2} x^2 + \frac{\theta}{2} \right) \psi \, ,\end{align} \tag{ 4.37 }$$

which we identify as the Euclidean Schrodinger equation $\partial_t \psi = - H_{\text{OU}} \psi$ where

$$\begin{align} H_{\text{OU}} = - \frac{\sigma^2}{2} \frac{\mathrm d^2}{\mathrm dx^2} + \frac{\theta^2}{2 \sigma^2} x^2 - \frac{\theta}{2} \, ,\end{align} \tag{ 4.38 }$$

in agreement with the result for H_OU we quoted earlier in (4.12). The operator (4.38) is now self-adjoint with respect to the inner product on twice-differentiable square-normalizable functions, as is appropriate for a Hamiltonian operator in quantum mechanics.

Combining equations (4.35), (4.36), and (4.38) gives the relation

$$\begin{align} \mathcal{L} = - \exp \left( - \frac{\theta}{2 \sigma^2} x^2 \right) H_{\text{OU}} \exp \left( \frac{\theta}{2 \sigma^2} x^2 \right) \, .\end{align} \tag{ 4.39 }$$

In terms of the stationary distribution (4.11) for the Ornstein–Uhlenbeck process, equation (4.39) can be expressed as

$$\begin{align} H_{\text{OU}} = - \mathbb{P}_s \left( x \right)^{-1/2} \, \mathcal{L} \, \mathbb{P}^{1/2}_s \,,\end{align} \tag{ 4.40 }$$

Thus, by performing the similarity transformation (4.40), we can convert from the (non-self-adjoint) Fokker–Planck operator $\mathcal{L}$ that generates time evolution of the probability distribution $\mathbb{P}(x, t)$ to the self-adjoint operator H_OU which generates time evolution of the wavefunction $\psi(x, t)$. Although we have reviewed the steps of this procedure for the Ornstein–Uhlenbeck process, the same argument applies to more general SPs which admit a stationary distribution, $H = - \mathbb{P}_s ( x )^{-1/2} \, \mathcal{L} \, \mathbb{P}^{1/2}_s$.

These observations suggest a numerical algorithm that one can perform to compute the energies and eigenstates associated with any SP that admits a stationary distribution. First, perform many draws from this stationary distribution to obtain an ensemble of initial positions x₀. Then use the update rule of the SP to time-evolve each initial position by a small time step Δt, which gives a collection of one-step sample paths $(x_0, x_{\Delta t})$. Bin these $(x_0, x_{\Delta t})$ pairs to create a histogram, which represents a discrete approximation to the joint probability distribution $\mathbb{P} ( ( x_{\Delta t}, \Delta t ) , ( x_0, 0 ) )$. We then normalize this histogram so that the sum of the entries in each row is unity; this yields a discretized estimate of the conditional probability distribution

$$\begin{align} \mathbb{P} \left( x_{\Delta t} , \Delta t \mid x_0 , 0 \right) = \frac{\mathbb{P} \left( \left( x_{\Delta t}, \Delta t \right) , \left( x_0, 0 \right) \right)}{\mathbb{P} \left( x_0, 0 \right) } \, ,\end{align} \tag{ 4.41 }$$

where the denominator $\mathbb{P} ( x_0, 0 )$ is the marginal distribution for the position x₀ at time 0. On the one hand, this conditional distribution obeys

$$\begin{align} \mathbb{P} \left( x_{\Delta t} , \Delta t \right) = \int\mathrm d x_0 \, \mathbb{P} \left( x_{\Delta t} , \Delta t \mid x_0 , 0 \right) \cdot \mathbb{P} \left( x_0 , 0 \right) \, ,\end{align} \tag{ 4.42 }$$

but on the other hand, as $\mathcal{L}$ generates infinitesimal time evolution of $\mathbb{P} ( x , t )$, we have

$$\begin{align} \mathbb{P} \left( x_{\Delta t} , \Delta t \right) \approx \mathbb{P} \left( x_0 , 0 \right) + \Delta t \, \mathcal{L} \, \mathbb{P} \left( x_0 , 0 \right) \, .\end{align} \tag{ 4.43 }$$

We now compare these two expressions to obtain an estimate for $\mathcal{L}$. In our discrete approximation, the integral in equation (4.42) becomes a matrix-vector product; the (i, j) entry of the matrix $\mathbb{P} ( x_{\Delta t} , \Delta t \mid x_0 , 0 )$ represents the conditional probability that $x_{\Delta t}$ lies in bin i given that x₀ lies in bin j. Therefore, in matrix notation, we have

$$\begin{align} \mathbb{P} \left( x_{\Delta t} , \Delta t \right)_i = \sum_j \mathbb{P} \left( x_{\Delta t} , \Delta t \mid x_0 , 0 \right)_{ij} \mathbb{P} \left( x_0 , 0 \right)_j \approx \mathbb{P} \left( x_0 , 0 \right)_i + \Delta t \, \sum_j \mathcal{L}_{ij} \, \mathbb{P} \left( x_0 , 0 \right)_j \, .\end{align} \tag{ 4.44 }$$

Solving for the matrix discretization of the Fokker–Planck operator $\mathcal{L}$ yields

$$\begin{align} \mathcal{L}_{ij} \approx \frac{\mathbb{P} \left( x_{\Delta t} , \Delta t \mid x_0 , 0 \right)_{ij} - \delta_{ij}}{\Delta t} \, ,\end{align} \tag{ 4.45 }$$

where the Kronecker delta δ_ij gives the entries in the identity matrix. This formula allows us to numerically compute the approximate Fokker–Planck operator. Next one performs the similarity transformation (4.40) using the stationary distribution; if $\mathbb{P}_s$ is not known in closed form, it can be estimated by simulating many sample paths with arbitrarily chosen initial positions, evolving each for a long time, and recording the steady-state distribution to which the ensemble settles down in the long run. Finally, one numerically diagonalizes the matrix H to obtain approximate energy eigenvalues and eigenfunctions.

The energy eigenvalues computed using this procedure are displayed in figure 4. Because the true spectrum of H_OU is $E_n = \theta n$, which is independent of σ, we only show results for different values of θ and fix σ = 1 across all experiments. The numerically estimated wavefunctions of the first three energy eigenstates are shown in figure 5 and compared to the true eigenstates. As we mentioned above, there are various sources of numerical error in this procedure, including Monte Carlo error, discretization and binning error which depend on Δt, and matrix diagonalization error in computing the eigenvalues of H; all of these sources contribute to the errors seen in figures 4 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Deep NN-QM with Ornstein–Uhlenbeck inputs

We now turn to an investigation of deep NN quantum mechanics, where a SP x_t is obtained via a NN whose inputs are themselves SPs satisfying all of the OS axioms. Since the Ornstein–Uhlenbeck process is the unique stationary Gaussian process that is also Markov, this is a natural place to begin, but we mention that one could also use other SPs as inputs. For instance, one could consider Feller processes like the Cox–Ingersoll–Ross model, defined by the SDE

$$\begin{align} \mathrm d x_t = - \theta x_t \, \mathrm dt + \sigma \sqrt{x_t } \,\mathrm d W_t \, ,\end{align} \tag{ 4.46 }$$

or the Jacobi diffusion,

$$\begin{align} \mathrm d x_t = - \theta x_t \, \mathrm dt + \sigma \sqrt{1 - x_t^2} \,\mathrm d W_t \, .\end{align} \tag{ 4.47 }$$

Both of these differ from the Ornstein–Uhlenbeck process in that the noise function which multiplies $\mathrm d W_t$ itself depends on x_t. Since the noise controls commutators and thus the effective $\hbar$ in the model, these processes can be roughly viewed as quantum systems with a ‘position-dependent Planck constant.’

We instead use standard OU processes $y_t^{(j)}$ and focus on the architecture

$$\begin{align} x_t = \frac{1}{\sqrt{N}} \sum_{i = 1}^{N} \sum_{j = 1}^{d} w_i^{\left(1\right)} \sigma \left( w_{ij}^{\left(0\right)} y_t^{\left(j\right)} \right) \, .\end{align} \tag{ 4.48 }$$

Like the choice of input SP, it would also be interesting to consider other architectures besides (4.48) and to investigate how architectural choices affect the physical properties of the quantum mechanical model.

In the numerical simulations below, for concreteness we take d = 3, N = 10, and we choose the parameters of the OU processes $y_t^{(j)}$ to be $\theta^{(j)}$ and $\sigma^{(j)}$. The noise parameters $\sigma^{(j)}$ are not to be confused with the activation function σ in equation (4.48), which we take to be

$$\begin{align} \sigma \left( z \right) = \tanh \left( z \right) \, .\end{align} \tag{ 4.49 }$$

One can either choose the weights $w_i^{(1)}$ and $w_{ij}^{(0)}$ to be identical across all sample paths x_t—in which case y_t is a deterministic function of a linear combination of RP processes, which is RP by propositions 2 and 3—or we can take the weights to be random variables which are drawn independently for each sample path, which is RP by the generalized result given in theorem 4.1. Here we take the latter approach; this is not a fundamental requirement of the deep NN-QM framework, and one could also study models with weights that are fixed across sample paths, but we choose the random-weight case since this illustrates the novel construction proposed in this work. Specifically, all of the weight variables are taken to be independent and identically drawn from unit Gaussians,

$$\begin{align} w \sim \mathcal{N} \left( 0 , 1 \right) \, ,\end{align} \tag{ 4.50 }$$

where w schematically represents all $w_i^{(1)}$ and $w_{ij}^{(1)}$.

Let us briefly comment on the bias variables, which we have set to zero in this construction. The reason is related to the mixing property, or cluster decomposition. As we alluded to above, if y_t is a mixing process and f is a measurable deterministic function, then $x_t = f ( y_t )$ is also mixing. However, if f is a random function, then mixing may not be preserved. A simple example is adding a random shift $\theta \sim P ( \theta )$ which is drawn separately for each sample path. If y_t is mixing and

$$\begin{align} \phi_\theta \left( y \right) = y + \theta \, ,\end{align} \tag{ 4.51 }$$

then the process $x_t = \phi_\theta ( y_t )$ is not mixing in general. An intuitive way to see this is by considering the centered covariance

$$\begin{align} \langle x_t x_s \rangle - \langle x_t \rangle \langle x_s \rangle & = \langle \left( y_t + \theta \right) \left( y_s + \theta \right) \rangle - \langle y_t + \theta \rangle \langle y_s + \theta \rangle \nonumber \\ & = \left( \langle y_t y_s \rangle - \langle y_t \rangle \langle y_s \rangle \right) + \left( \langle \theta^2 \rangle - \langle \theta \rangle^2 \right) \, ,\end{align} \tag{ 4.52 }$$

where we have used that θ and y are uncorrelated to cancel the cross-terms. The quantity in the first set of parentheses is the centered covariance of process y_t while the second set of parentheses contains the variance of θ. Thus, even if the centered covariance of y_t decays to zero for large values of $| t - s |$, the corresponding covariance of process x_t approaches the constant value $\mathrm{var} ( b )$, which is non-zero if b is a random variable.

Such a random shift (4.51) is essentially what randomly chosen bias variables in the NN implement, and indeed we see that the mixing property (or cluster decomposition) fails if we include these random biases in our network. Therefore, we turn off the bias parameters, which leads to an output process y_t whose empirical two-point function obeys cluster decomposition, as we investigate below.

Two-point function

As in our discussion of the Ornstein–Uhlenbeck process, the first natural observable to study for our deep NN-QM model is the two-point function. Because we take d = 3 input OU processes $y_t^{(i)}$, each of which is determined by a drift parameter $\theta^{(i)}$ and a noise amplitude $\sigma^{(i)}$, there are six available parameters to tune in the definition of x_t. The experimental values of $G^{(2)} ( t, 0 )$ for several choices of these parameters are displayed in figure 6. All other parameters in the NN—including the width N = 10 of the hidden layers, the distribution (4.50) of weights, and the choice (4.49) of tanh activation function—are fixed across all of the experiments.

Zoom In Zoom Out Reset image size

The results of these simulations confirm several basic properties that one might expect of our deep NN-QM model. First, the correlation function $G^{(2)} ( t, 0 )$ decays to zero at large t, which gives evidence that mixing or cluster decomposition holds in this theory (as we mentioned above, this is a consequence of our choice to set the biases to zero, and mixing fails if we choose random non-zero biases). Second, we see that the behavior of the two-point function changes as the parameters $\theta^{(i)}$ and $\sigma^{(i)}$ of the input OU processes are varied, as it should. Third, although we have no theoretical prediction for an analytical expression for the two-point function in this model to which the experimental results could be compared, by equation (2.5) we expect the leading exponential fall-off of $G^{(2)}$ to be determined by the gap E₁ from the ground state to the first excited state (since E₀ = 0 by construction in the stationary state of a SP). By performing exponential fits with the functional form $G^{(2)} = A \mathrm e^{- \lambda t}$, we therefore obtain approximations of the first excited state energies in these deep NN-QM models, which are also recorded in figure 6.

Commutator

For the Ornstein–Uhlenbeck process, we found that the discrete approximation C(t) of the commutator $[ \hat{x} ( t ) , \hat{p} ( t ) ]$, defined in equation (4.20), was constant in time with a magnitude set by the effective $\hbar$ of the SP. This should, of course, hold for any conventional quantum system with a quadratic kinetic term. However, since the deep NN-QM that we have defined here is obtained from a fairly involved transformation of several OU processes, it is not clear a priori that it corresponds to a familiar quantum system for which $p(t) \sim \dot{x} ( t )$ so that C(t) is constant.

Nonetheless, repeating the numerical computation of C(t) via simulations as described above, we find that this is indeed the case. These results are presented in figure 7, which shows that C(t) is also constant for deep NN-QM, with an effective $\hbar$ that is a function of the noise parameters $\sigma^{(i)}$ defining the input OU processes.

Zoom In Zoom Out Reset image size

We note that the crucial non-differentiability of paths, which is required for non-trivial commutation relations, in the deep NN-QM approach is inherited from the non-differentiability of the input processes. However, the character of this non-differentiability is ‘re-shaped’ by the application of the NN. One way to see this is to note that the commutator values in the transformed process of figure 7 differ from those in the input processes, even in the case where all σ_i are equal. Another simple way to see this is to consider a deterministic NN $\phi_\theta ( y_t ) = y_t^2$. The SP obtained by squaring the Ornstein–Uhlenbeck process is related to the Cox–Ingersoll–Ross model (4.46), whose noise (and hence non-differentiability) is position-dependent, so even in this toy example qualitative aspects of the roughness of paths are considerably modified.

Uncertainty relations

It is also interesting to consider the behavior of the uncertainties Δx and Δp in the deep NN-QM model. These uncertainties must, of course, obey the Heisenberg uncertainty relation, but we no longer expect the product $\Delta x \, \Delta p$ to saturate the bound, since this is a special property of minimum uncertainty states such as Gaussian wavepackets.

The simulated uncertaintes $\left( \Delta x \right)^2$ and $\left( \Delta p \right)^2$ for several choices of parameters are displayed in table 2. Again, we have no theoretical prediction for the product of uncertainties, but the lower bound implied by the Heisenberg uncertainty relation is displayed in the final column of table 2. This lower bound is not uniform across all rows of this table since it is set by the effective value of $\hbar$ in the model. We estimate $\hbar_{\text{eff}}$ by computing the commutator C(t), and as we have already seen in figure 7, this quantity depends on the parameters. We find that the Heisenberg uncertainty relation is not saturated for any choices of parameters that we consider in this deep NN-QM theory.

Table 2. We compute the uncertainties in position and momentum for the deep NN-QM process with three OU inputs. The definitions (4.28) and (4.25) are used to compute the momentum uncertainties. We simulate 200 epochs with $10\,000$ paths per epoch, each with 1000 update steps, for each combination of parameters, and the standard deviation across epochs appears following the symbol ± for each quantity. We compare the simulated uncertainties for several values of the parameters $\theta^{(i)}$ and $\sigma^{(i)}$ defining the input OU processes, and for two values of the time step $\mathrm dt$. We also include the value $\frac{\hbar_{\text{eff}}^2}{4}$, which gives a lower bound for the product $\left( \Delta x \right)^2 \cdot \left( \Delta p \right)^2$ by the Heisenberg uncertainty principle, where $\hbar_{\text{eff}}$ is estimated from the commutator C(t) which approximates $[\hat{x} , \hat{p} ]$. We find that $\left( \Delta x \right)^2 \cdot \left( \Delta p \right)^2 \gt \frac{\hbar_{\text{eff}}^2}{4}$ in all cases, so the uncertainty bound is not saturated.

θ1	θ2	θ3	σ1	σ2	σ3	$\mathrm dt$	$\left( \Delta x \right)^2$	$\left( \Delta p \right)^2$	$\left( \Delta x \right)^2 \cdot \left( \Delta p \right)^2$	$\frac{\hbar_{\text{eff}}^2}{4}$
$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	0.100	0.3031 ± 0.0029	0.1360 ± 0.0018	0.0412 ± 0.0008	0.0307 ± 0.0005
$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	0.010	0.2997 ± 0.0033	0.1480 ± 0.0178	0.0444 ± 0.0053	0.0307 ± 0.0006
1	1	1	1	1	1	0.100	0.4175 ± 0.0037	0.9227 ± 0.0086	0.3852 ± 0.0066	0.2687 ± 0.0044
1	1	1	1	1	1	0.010	0.4104 ± 0.0040	1.2741 ± 0.0553	0.5229 ± 0.0232	0.2850 ± 0.0052
2	2	2	2	2	2	0.100	0.5393 ± 0.0044	4.8038 ± 0.0352	2.5908 ± 0.0387	1.9196 ± 0.0292
2	2	2	2	2	2	0.010	0.5236 ± 0.0047	10.4247 ± 0.1741	5.4586 ± 0.1127	2.4223 ± 0.0360
$\frac{1}{2}$	1	2	$\frac{1}{2}$	1	2	0.100	0.4365 ± 0.0042	2.1579 ± 0.0235	0.9420 ± 0.0183	0.6104 ± 0.0121
$\frac{1}{2}$	1	2	$\frac{1}{2}$	1	2	0.010	0.4263 ± 0.0039	3.7006 ± 0.1075	1.5777 ± 0.0515	0.6834 ± 0.0150
2	1	$\frac{1}{2}$	2	1	$\frac{1}{2}$	0.100	0.4371 ± 0.0046	2.1610 ± 0.0255	0.9446 ± 0.0204	0.6100 ± 0.0120
2	1	$\frac{1}{2}$	2	1	$\frac{1}{2}$	0.010	0.4255 ± 0.0039	3.6972 ± 0.1051	1.5731 ± 0.0492	0.6819 ± 0.0138
$\frac{1}{2}$	2	1	2	$\frac{1}{2}$	1	0.100	0.5329 ± 0.0059	1.0278 ± 0.0113	0.5477 ± 0.0108	0.2674 ± 0.0054
$\frac{1}{2}$	2	1	2	$\frac{1}{2}$	1	0.010	0.5287 ± 0.0061	1.8866 ± 0.0762	0.9974 ± 0.0415	0.3101 ± 0.0066
1	$\frac{1}{2}$	2	1	2	$\frac{1}{2}$	0.100	0.5330 ± 0.0056	1.0282 ± 0.0107	0.5481 ± 0.0105	0.2671 ± 0.0053
1	$\frac{1}{2}$	2	1	2	$\frac{1}{2}$	0.010	0.5286 ± 0.0057	1.8828 ± 0.0840	0.9952 ± 0.0451	0.3098 ± 0.0057
2	$\frac{1}{2}$	2	$\frac{1}{2}$	1	1	0.100	0.3678 ± 0.0038	0.9056 ± 0.0099	0.3331 ± 0.0063	0.1943 ± 0.0035
2	$\frac{1}{2}$	2	$\frac{1}{2}$	1	1	0.010	0.3603 ± 0.0042	1.1025 ± 0.0492	0.3973 ± 0.0188	0.1934 ± 0.0037

Spectrum and eigenfunctions

Finally, we repeat the calculation of the energy eigenvalues and eigenfunctions for the deep NN-QM model. These results are shown in figure 8. We see that the energy eigenvalues are approximately linear in n for small n, as in the harmonic oscillator, but upon investigating higher energy eigenvalues we find that the curve is concave down. The rate at which the energies ‘level off’ seems to be set by the θ parameters of the input Ornstein–Uhlenbeck processes, with larger θ values giving models that deviate from the linear SHO-like spectrum more quickly. We also see that, unlike the OU example, the spectrum appears to depend on the noise parameters $\sigma^{(i)}$ in addition to the $\theta^{(i)}$.

Zoom In Zoom Out Reset image size

We do not display the energy eigenfunctions since their shape is very similar to that of the OU eigenfunctions in figure 5. This suggests that—at least for low-lying states—the effect of the deep NN transformation primarily modifies the energy eigenvalues, while leaving the energy eigenfunctions essentially unchanged. This is similar to the behavior of deformations of quantum systems where the seed Hamiltonian H₀ is mapped to a function $H = f ( H_0 )$ of the undeformed Hamiltonian, such as in the quantum-mechanical versions of the ${T\bar{T}}$ [43, 44] and root-${T\bar{T}}$ [45, 46] deformations. However, at high energies, we see that the deep NN-QM spectrum levels off, so the effect of the transformation on the eigenfunctions for high-energy states may be more dramatic.

Although we will not undertake a detailed sensitivity analysis for deep NN-QM in this work, we note that for each of the quantities considered here, one could also vary other pieces of defining data in a deep NN-QM model and study how the resulting observable changes. For instance, one could promote the commutator C(t) displayed in figure 7 to a family $C ( N, t )$ of commutators, where N is the number of neurons in the architecture (4.48), and study the sensitivity of this quantity to N. Since modifying the architecture changes the theory defined by the deep NN-QM, such an analysis could be interpreted as studying how observables vary in families of quantum mechanical models.