Quantum Hamilton equations for multidimensional systems

Within a stochastic picture quantum systems can be described in terms of kinematic and dynamic equations where the particles follow a conservative diffusion process. Pavon has shown that these equations and the Schrödinger equation can be derived from a quantum Hamilton principle, generalizing the Hamilton principle from classical mechanics. As shown recently by Köppe et al for the one-dimensional case, the reformulation of the quantum Hamilton principle as a stochastic optimal control problem leads to quantum Hamilton equations of motion which can be seen as a non relativistic generalization of Hamilton’s equations of motion to the quantum world. In this report we will extend this to higher dimensional problems and derive, by incorporating the concept of supersymmetric quantum mechanics, a general scheme that allows the determination of all bound states within the stochastic framework. Physical principles from classical mechanics like the decoupling of the center-of-mass motion in multi-particle systems or the solution of the Kepler problem, as a special case of the two-body problem, using space-time symmetries, are translated to this formulation of quantum mechanics. We will present numerical results for the ground state and lower excited states of the hydrogen atom in Cartesian coordinates, and, additionally, we will show that for spherically symmetric potentials the solution of the two-body problem can be reduced to a system of equations concerning only the radial part.

Keywords: quantum mechanics, stochastic mechanics, hydrogen atom (Some figures may appear in colour only in the online journal)

Introduction
Unlike any other theory in physics, quantum mechanics gave and still gives a lot of room for discussion. This is due to the fact that it conceptually differs from classical mechanics in its standard formulation and that we are concerned with phenomena that can not be described based on everyday experience any longer. This resulted in the need for a new physical theory and with that two main components that had to be put together: (1) the mathematical formalism and (2) its physical interpretation. The most famous and accepted formalism relies on the Schrödinger picture where the time evolution of a quantum system is characterized by a complex wave function ψ in Hilbert space, which solves a linear partial differential equation, the Schrödinger equation [1]. However, since ψ is an abstract quantity, its physical interpretation becomes involved because it is not as natural or intuitive as the objects of classical mechanics. For example, its meaning only became clear when Born postulated that |ψ(x, t)| 2 is the probability density of observable x. This is one of the reasons why there are so many different approaches to describe quantum systems. In 1966, Edward Nelson added his stochastic picture of quantum mechanics [2] to the existing pool of approaches. He was able to establish a set of forward-backward stochastic differential equations (FBSDEs) as a time-reversible description of the stochastic path of a particle. The Schrödinger equation is then derived requiring the mean acceleration to be equal to the negative gradient of the potential divided by the mass m of the particle, and by stating that the diffusion constant, denoted by σ 2 , equals m . Within this stochastic setup, a lot of aspects come naturally, e.g. that the probability to find the particle at position x at time t is given by ρ(x, t)dx = |ψ(x, t)| 2 dx is intrinsic to this formalism, where ρ is the probability distribution associated with the stochastic process of the particles path. The stochastic mechanics approach allows to study quantum mechanical systems from a different point of view and can give seminal results [3][4][5]. However, the problem within Nelson's stochastic picture is that it is inconvenient, or, in general even not possible, to find solutions to his FBSDEs, and thus one was limited to solving the Schrödinger equation beforehand.
In [6] it was possible to derive modified FBSDEs for unidimensional, non relativistic, non stationary quantum systems with the help of stochastic optimal control theory and in [7] it was shown that the combination of this approach with tools from the concept of supersymmetric quantum mechanics (SUSY) allows for a determination of the complete bound spectrum (energies and eigenstates) of 1d quantum systems. The intention of this article is to show that we can generalize this formalism to multiple dimensions and analyse quantum problems using symmetries and conservation laws as in classical mechanics. This is exemplified by solving the two-body problem (hydrogen atom) with the quantum Hamilton equations. Additionally, this formalism is not meant to supplant or outdo the standard techniques used to solve the Schrödinger equation. The aim rather is to extend the number of methods available for and the scope of the treatment of quantum mechanical problems in general.
The following section is devoted to the derivation of the quantum Hamilton equations. In the third section, the fact that one deals again with particles in a classical sense will be exploited by showing that one can follow the same physical arguments in the solution of the quantum problem as in classical mechanics. This will be shown concerning the decoupling of the center-of-mass (COM) motion of multi-particle systems. The last section will deal with the two-body problem, namely the Kepler problem. A part of this is dedicated to the incorporation of the SUSY approach into the stochastic formulation, which allows us to determine the excited states of a quantum system iteratively from the ground state solution of the quantum Hamilton equations. We will provide numerical results for the two-body problem and propose a way to reduce problems of spherical symmetry to the solution of the radial part only. In particular, this will be discussed for the ground and excited states of hydrogen atomlike systems.

Quantum Hamilton equations
Assume that we are concerned with a system of N particles moving in d dimensions characterized by their position x(t) at time t ∈ [0, T] where (x t ) = (x(t)) 0 t T is a stochastic process in R N·d=:m . In that case, the quantum Hamilton principle according to Pavon [8] can be written as a stochastic optimal control problem where the cost functional is to be extremized w.r.t. the quantum velocity R m v q = v − iu. Here E denotes the expectation value w.r.t. the probability density ρ(t, x(t)) associated to (x t ), L = T − V is the (stochastic) Lagrangian with kinetic energy T and potential energy V of the system and the terminal value The quantum velocity is the optimal feedback control to x(t) and the cost functional is subject to the controlled (doubly) stochastic differential equation (2) This equation can be derived by combining the FBSDEs of Nelson [2]. Here σ is a tensor depending on the masses of the N particles W(t) 0 t T is a m dimensional forward in time Wiener process and W − (t) 0 t T is a m dimensional backward in time Wiener process. The cost functional in equation (1) is a combination of two saddle-point variational principles in its real and imaginary parts, respectively, referring to the action and the entropy production of the system [8]. In order to save space, short hand notations like v q (t, x) = v q (t, x(t)) were and will be used throughout this report.
For the control problem in [8] we can construct the associated stochastic Hamiltonian [9,10] Analogously to the deterministic maximum principle, adjoint variables λ(t) ∈ C m and q(t) ∈ C m×m occur. This is due to the stochasticity of the terms in the expectation value of the cost functional, where one has to introduce random variables that find the optimal trade-off concerning the uncertainty of the system. These variables satisfy the following BSDEs [9,10] Finding the roots of the Hamiltonian w.r.t. v q relates the adjoint processes with the quantum velocity where m represents the mass tensor of the system and ∂ x v q denotes the Jacobi matrix of v q w.r.t. x. The last equation results from Itô's formula applied to v q (x, t). By separating the real and imaginary parts of those equations and defining the (real) physical one ends up with the following system of coupled FBSDEs which are called the quantum Hamilton equations. Note that the first two equations are Nelson's FBSDEs that were used to derive (2). The sanity check of those equations by taking the classical limit (σ → 0, |u| → 0) leads directly to Hamilton's deterministic equations of motion. The structure is the same as in classical mechanics which enables us to use similar logic and approaches for quantum systems within this framework. This will be discussed in the following section by having a look at multi particle systems.

Quantum Hamilton equations for multi particle systems
In the previous section we saw that there exists a kinematic and dynamic mathematical description of non relativistic quantum systems as a generalization of classical mechanics. This leads immediately to the question whether we can draw some analogies, regarding properties and also concepts, to classical analytical mechanics. We take as an example the QHEs for many particle systems without external potential, leading to similar properties as in the classical description. For the sake of simplicity, let us consider the stationary case and bound solutions first. This results for ground state wave functions in v = 0 and u = u(x(t)) [11]. However, it is also possible to consider the following for non-stationary systems as well. We are concerned with the special case of an atom-or ion-like system where we deal with a nucleus of mass m 0 and n electrons with masses m denoted by R 0 ∈ R d and R i ∈ R d (i ∈ {1, . . . , n}), respectively. The cost function in the stationary case reads [6] and it is subject to the kinematic equation is the optimal feedback control for particle i and σ i = mi .
Note that the backward SDEs are omitted here, but they are still (additional) constraints.
When thinking of multi-particle systems in classical mechanics, the first thing that comes to mind is that the total momentum is constant due to spatial translation invariance. Hence, we introduce the relative coordinates of every electron w.r.t. the nucleus and the COM, respectively, Here M = m 0 + n m stands for the total mass of all particles. The osmotic velocities should be defined accordingly, as (11) The above defined relative coordinates and COM coordinates are again stochastic processes. They obey where the drift terms and diffusion constants can be obtained using Itôs formula. With the new stochastic processes R and r i and the definition of the osmotic velocities of the COM U and the relative coordinates u i the maximum cost functional, particularly the Lagrangian in (8) has to be adjusted to the new coordinates. This yields the Lagrangian where µ := m M (M − m). By taking into account that the potential V is an interaction potential which only depends on distances between the particles, i.e. r 1 , r 2 , . . . and r ij = |r i − r j |, equation (14) separates the part depending on U 2 from the one containing the u 2 i . If we define stochastic processes and matricesm the amount of terms of the associated Hamiltonian is condensed. Omitting the time dependence, e.g. replacing X = X(t), we end up with At this point point we have to introduce three different stochastic processes (p(t)), (λ(t)) ⊂ R (n+1)·d and (q(t)) ⊂ R (n+1)·d×(n+1)·d in appropriate dimensions with p(t) as a backward in time process adapted to the FSDE of (X(t)) and (λ(t)) as a forward in time process adapted to the BSDE of (X(t)). The processes λ and p satisfy certain SDEs [10,12] and can be connected to the feedback control Y by extremizing the Hamiltonian H with respect to Y. This results in SDEs for a new process Ȳ as follows Recalling that the Potential does not depend on R the SDEs for R can be written as The second equation has only terms including Wiener processes. If we write q again in terms of Jacobi matrices with respect to x, we can conclude that U ≡ U T = const is a solution to the problem. If we take the expectation value of (16) and of the backward Recalling that we're interested in bound solutions to this problem, formally this means that there is no stationary state for the center of mass motion. If we consider the non-stationary case, the presented transformation would be quite similar with analogous definitions for the current velocity of the COM, V , and the relative coordinates, v. Then we would also have U ≡ 0, but V ≡ V 0 = const. This refers to a free particle in the quantum case, e.g. the plane wave ψ ∝ exp[ i M V 0 · R] and in the classical limit to the conservation of the total momentum. Last but not least, it is also possible to state that the expectation value of the total angular momentum defined as . The calculation is similar to classical physics with the difference that dL total leads to terms where only products including dW (−) occur.
For the relative motion, the QHEs are obtained as follows (22) leave us with a system of 2 · d · n coupled SDEs where the COM motion is separated. This is in analogy to classical mechanics, where the conservation of the total momentum allows to separate this motion from the relative coordinates. In addition, this shows also that analytical or even numerical solutions within the stochastic approach are limited to few-body problems highlighting the close relationship of the QHEs to the classical ones. In classical mechanics, one proceeds to solve the two-body central force problem exactly, i.e. the Kepler problem. In quantum mechanics this is the exactly solvable hydrogen atom. This will be part of the next section.

Hydrogen atom (two-body problem)
The hydrogen atom defines a bound state of a proton and electron with masses m n and m e and opposite charges e and −e (e > 0). In the non relativistic case, the interaction between the two particles is described by the Coulomb potential where e 2 0 = e 2 4πε0 and ε 0 is the permittivity of vacuum. In the case of a two-body problem, the previous section leaves us with dynamical equations for R leading to a (free) 3d Brownian motion of the center of mass. This is different for r. According to (22), the SDEs for the relative coordinate read This is a Kepler-like system for the relative coordinate. This becomes more clear when we consider the expectation value which corresponds to Ehrenfest's theorem describing the time derivative of the expectation values of the position and momentum operators in relation to the gradient of the potential. Since we already know from the ground state solution of the Schrödinger equation that there is no angular momentum in this state, it is interesting to ask why the particle in the stochastic picture does not get trapped by the proton. This would be the case for the classical Keplerproblem since there some angular momentum would be needed that prevents the particle from falling into the nucleus. It is possible to solve these SDEs for the relative coordinate via direct evaluation of the SDEs forward and backward in time by discretizing time and spatial axis and a least squares method for conditional expectations [6].
The advantage of that description is that after determining the velocity fields, in this case only the osmotic velocity, we can visualize possible trajectories that a particle may follow.

Excited states
Recently, we showed that it is possible to construct excited states for stationary systems in the frame of Nelson's stochastic mechanics in one dimension using the concept of supersymmetry and the quantum Hamilton equations [7]. The purpose of this section is to show how this can be generalized to multi-dimensional and multi-particle systems.
The concept of supersymmetry in one-dimensional quantum mechanics introduced in 1982 by Witten [13] is a well-known concept in mathematics. In 1912, Darboux [14] derived a method for Sturm-Liouville type differential equations to construct an infinite series of exactly solvable problems starting from a known base problem in one dimension. In 1955, Crum [15] presented a compact form of the expressions and investigated the Schrödinger equation. The work of Sabatier [16] generalizes this idea to multi-dimensional systems and a lot of work in multi-dimensional SUSY quantum mechanics followed [17].
In our case, the strategy is to transform the well-known formalism from wave functions to the velocity fields of Nelsons equation by starting with the stationary case where the quantum mechanical Hamiltonian reads with p 0 ( p 0 1) bound states and x = (x 1 , ..., x d ) T ∈ R d . Note that atomic units are used for simplicity and that ∂ i ≡ ∂ xi . To simplify the calculations, the potential is shifted by V 0 (x) → V 0 (x) + const to have zero ground state energy. In order to proceed we will state the following: if H 0 has a normalized ground state wave function, then, in any dimension, this function has zero nodes [18]. It's possible to find a well defined function χ(x) to write the ground state wave function in the following form The superscript (n) in ϕ (n) k indicates the nth excited state of the wave function belonging to the Hamiltonian H k . Note that n has to replaced by a vector in multiple dimensions for excited states.
To factorize the Hamiltonian H 0 define the 2d operators fulfilling Consequently, the Hamiltonian can be factorized to (using Einstein's sum convention) If we want to apply the Darboux transformation to the Schrödinger equation with Hamilton operator H 0 , the partial differential equation can be written in terms of a partner Hamiltonian (all indices take d values) Despite the fact that the result is a Hamiltonian in matrix form and it could be a difficult task to find eigenfunctions, there exist some practical intertwining relations. Start with The diagonal elements (i = k) fulfill Following the procedure of [7] the relation is the (n + 1)th excited wave function in the quant um number corresponding to the kth coordinate and ϕ (n),k 1 is the nth state for H k,k 1 . The ground state of H k,k 1 is given by the QHEs (7) with a modified potential V k,k where in the stationary case for ground states we have v = 0 and q 0 (t) = σ∂ x u 0 (x). The calcul ation of the excited states is reduced to the determination of ground states of the partner Hamiltonians. Due to the node theorem these states are node-free and therefore their corresponding velocities have no singularities which in turn leads to diffusions covering the whole space. That allows to determine the original excited wave functions with their nodes by the application of differential operators. In the end this leads in an iterative way to all bound states ϕ (a1,...,b k ,...,c d ) 0 for H 0 . For example, to calculate the state ϕ . From here we can determine a first excited state velocity w.r.t. the axis x 1 , e.g. and finally, u (2),(1,1,0) 0 similar to (41). Thus, this iterative method allows us to calculate the excited states by simply adjusting the potential with the help of previously calculated velocities and, finally, by acting with operators on the newly determined osmotic velocities. The solutions of the QHEs in (40) can be found according to a numerical iteration similar to the presented algorithm in [6] but now in the case of a multi dimensional problem.
As an example we take again the two-body problem from the previous section and determine one first excited state in Cartesian coordinates. Once more, two paths of a first excited state are shown on the rhs of figure 1, where the motion hints at the typical dumbbell shaped p -orbital. The numerical expectation value in this case gives for the energy (in Rydberg units) E 1 = −0.12 ± 0.08 and for the angular momentum it gives L x ϕ (1,0,0) 0 = 0.98 ± 0.08, where the errors are, again, given w.r.t. five independent simulations. Thus, the first excited state in x direction can be determined to be the p x orbital.
We see that the numerical solutions give good estimates in comparison to the exactly known solution. However, the problem with the excited states (and also the ground state) is, that a lot of numerical effort is needed due to a 2 · 3 dimensional problem that is to be solved within a numerical approach including numerical derivatives. So we run into the typical problems of many particle physics whether it is classical physics where the size of phase space explodes or quantum mechanics where the configuration space explodes with increasing number of particles. That is the reason why there is a need of a new algorithm or different approach if we want to solve these kind of equations for higher dimensionality within that formalism. However, in the case of the hydrogen atom, the considered problem can be further simplified.

Using the symmetry of the hydrogen atom
If we stay with the hydrogen atom, another possible question that arises is whether it is possible to derive kinematic and dynamic equations within the framework of Nelson's stochastic mechanics for the ground state of the two-body problem where the problem should be reduced to a 2 · 1 dimensional problem where the radial part is separated from the angular part. Since the two-body problem possesses a spherical symmetry we may introduce spherical coordinates (r, ϑ, ϕ) ∈ (0, ∞) × (− π 2 , π 2 ) × [0, 2π) =: O. Using (24) the new coordinates r(r(t)), ϑ(r(t)) and ϕ(r(t)) obey the following the SDEs (using Itô's formula) The transformed three dimensional Wiener process W = SW is also a 3 dimensional Wiener process where S is the spherical transformation matrix. Ignoring the part concerning the COM motion, the Lagrangian reads In order to set up the stochastic Hamiltonian H r (r, ϑ, ϕ, u r , u ϑ , u ϕ ) we definẽ and get again with the three stochastic processes (p(t)), (λ(t)) ⊂ R 3 (backward/forward processes adjoint to the drift in the FSDE/BSDE) and (q(t)) ⊂ R 3×3 . It turns out, that λ(t) can chosen to be 0 ∀t ∈ [t 0 , T] and that the remaining stochastic processes obey dp = −∂ r Hdt + qdW − where, by definition, we here have ∂ r =r∂ r +θ∂ ϑ +φ∂ ϕ . The components of the process p = ( p r , p ϑ , p ϕ ) T are connected to the osmotic velocity by ( p r , p ϑ , p ϕ ) = µ(u r , u ϑ , u ϕ ). This results from maximizing H with respect to u. In the end this yields with the matrix valued stochastic process q =   q rr q rϑ q rϕ q ϑr q ϑϑ q ϑϕ q ϕr q ϕϑ q ϕϕ   and qdW − = There are no potential terms in the latter two SDEs. Now assume that u ϑ (r(t), ϑ(t), ϕ(t)) = 0 and u ϕ (r(t), ϑ(t), ϕ(t)) = 0 ∀t ∈ [0, T) with final values u ϑ (r(T), ϑ(T), ϕ(T)) = u ϕ (r(T), ϑ(T), ϕ(T)) = 0. If we recall that the matrix valued stochastic process q has to ensure that the pair (p(·), q(·)) or in this case (u(·), q(·)) is a solution to (49) it follows immediately that q ij (t) vanishes for i ∈ {ϑ, ϕ}, j ∈ {r, ϑ, ϕ}. Hence, if we insert this into the equation (49) this leaves us with a set of equations to solve for the radial part only, Thus, we derived QHEs that allow to solve a reduced one dimensional problem concerning the distance r between the proton and the electron for the ground state of the hydrogen atom. Note that this derivation is valid for any central symmetric potential, not only the Coulomb potential. The differences to the 1D problems so far are the additional drift terms which occur due to the transformation of variables. For example the term µr in (51) is a 'probabilistic' force which tries to push the particle away from the center because the volume element dV = r 2 sin ϑdrdϑdϕ gets smaller when r goes to zero. In the case of the Schrödinger equation treatment for spherical problems, these additional terms occur as consequence of the transformation of the Laplacian, whereas in the stochastic formulation they follow from the non vanishing variance of the stochastic process concerning the position in the mean square limit. This is also the answer to the question posed earlier why the electron does not simply fall into the proton. A numerical solution for (51) and (52) can be found using the algorithm mentioned in [6]. The result for the ground state of the hydrogen atom is shown in figure 3 (black). We see that the numerical result in that case is in agreement with the exact analytic solution for the ground state. Note that the motion of the electron in the ground state is diffusive while being tied to the proton as shown in the left plot of figure 1. At first this may seem odd, since one could expect a perturbed Kepler-like elliptic motion of the electron or something that may be related to the Bohr model of the atom. In the ground state, however, this is not observed due to the lack of an angular velocity when taking expectation values of the angular velocities. The classical Kepler-problem should be recovered in the Bohr correspondence limit which was shown for the Nelson diffusion at least in two dimensions [19].
As already shown in figure 3, the excited states can be determined as well, and in this case computationally easier than in Cartesian coordinates. The equations proposed in the previous section have to be adjusted since we consider the radial part of the problem in spherical coordinates. However, if we define a new osmotic velocitỹ which is exactly the drift term in (51), we can use the scheme as suggested before. This is in analogy to the transformation of the radial wave function R (r) = rR(r) when solving the Schrödinger equation that allows to find a (SUSY-)decomposition of the adjusted Hamilton operator. Note that in this case u r = µ ∂ r ln R. Let us denote the potential of the considered problem as V 0 , which is the Coulomb-potential. Then the first partner potential reads where u (0) 0,r is the radial osmotic velocity of the ground state concerning the potential V 0 . If we use the exact solution u (0) 0,r = − µa0 with a 0 as the Bohr radius, V 1 is exactly the effective potential with angular momentum quantum number l = 1. Recall that solving the Schrödinger equation and separating the radial part form the angular parts leads to the effective potential for the radial part. This means that the first partner potential gives the spectrum for l = 1. It should be noted that the ground state energy of this potential is equal to the first excited state of V 0 . This is also valid for the second partner potential r,0 of the hydrogen atom (n = 1, l = 0), the QHEs with the partner potential V 1 give the non-singular velocities of a node-free ground state u (1) r,0 concerning l = 1 states (n = 2). Using these two determined velocities the first excited state u (0) r,1 of V 0 is calculated (n = 2, l = 0). This state is depicted (orange) in figure 3. The iterative determination of the remaining excited states is then done by analogy. In general, it is possible to state that the ground state velocities for the partner potential V l = V 0 + l(l + 1) 2 2µr 2 (l ∈ N 0 ) read u (l) r,0 = − (l+1)µa0 + l µr . These ground-state velocities are non-singular which corresponds to node-free (radial) wave-functions. The energy of these states can be calculated via the expectation value w.r.t. t of the radial Hamiltonian as follows

Conclusion
We presented the derivation of quantum Hamilton equations for N particles in arbitrary dimensions extending the formalism presented in [6]. Using the SUSY approach to quantum mechanical problems, it is in principle possible to determine the entire bound spectrum of solutions. These quantum Hamilton equations allow a more intuitive description and physical interpretation of quantum systems because they are closely related to classical mechanics. Due to this analogy, one encounters similarities to classical problems in the stochastic picture. Quantum mechanical conservation laws of total momentum and angular momentum follow from the invariance of the QHEs under spatial translation or rotation in the same way as in classical mechanics. In particular, it was possible to reduce spherically symmetric problems within the suggested stochastic optimal control problem to a one dimensional radial problem. The algorithm proposed in [6] was applied successfully to the hydrogen atom for this reduced case. Furthermore, we set up a SUSY scheme for the hydrogen atom, which allows for the calculation of the radial part of all bound states and their energies.