Bouncing oil droplets, de Broglie's quantum thermostat and convergence to equilibrium

Recently, the properties of bouncing oil droplets, also known as"walkers", have attracted much attention because they are thought to offer a gateway to a better understanding of quantum behaviour. They constitute indeed a macroscopic realization of wave-particle duality, in the sense that their trajectories are guided by a self-generated surrounding wave. The aim of this paper is to develop a phenomenological theory for the behavior of walkers in terms of de Broglie-Bohm and Nelson dynamics. We study in particular how modifications of the de Broglie pilot-wave theory, \`a la Nelson, affect the process of relaxation to quantum equilibrium, and prove an H-theorem for the relaxation to quantum equilibrium under Nelson dynamics. We compare the onset of equilibrium in the Nelson and de Broglie-Bohm approaches and we also propose some simple experiments by which one can test the applicability of our theory to the context of bouncing oil droplets.


Introduction
"Walkers" are realized as oil droplets generated at the surface of a vibrating oil bath. As shown by Couder and Fort [1][2][3], the vibration of the bath prevents the coalescence of the droplets with the surface, allowing them to remain stable for very long times. Moreover, the trajectories of the walkers are guided by an external wave [4,5] that they themselves generate at the surface of the oil bath. From this point of view, walkers are reminiscent of wave-particle duality [2,6], and they seem to offer deep analogies with de Broglie-Bohm particles [7]. Up until now, different aspects of walker dynamics have been studied in a purely classical framework, typically in a hydrodynamical approach [3,5]. For instance, certain models address their deformations due to their bouncing off the surface of the bath, in function of the density and viscosity of the oil and other parameters [5]. Other studies describe the dynamics of the surface waves that the walkers generate during the bouncing process, and how those waves in turn guide their trajectories. In these models, this complex behavior is characterized by a memory time which relates the dynamics of the walker bouncing at time t to its successive bouncing positions in the past [8,9]. The presence of such a memory effect establishes a first difference with quantum mechanics. Normally, in quantum mechanics, it is assumed that all results of any possible future measurements to be performed on a quantum system are encapsulated in its present quantum state [10]: its wave function at the present time t.
Droplets also transcend the most common interpretations of quantum theory which prohibit any description of the system in terms of instantaneous, classical-like trajectories. Droplets and their trajectories are visible with the naked eye at any time and standard interpretations of quantum mechanics do not apply. This is why we believe that it is necessary and worthwhile to adapt realist (causal) formalisms such as de Broglie-Bohm (dBB) dynamics [11,12] or a stochastic version thereof à la Nelson [13] to explore the analogy with quantum systems. This is the main motivation of the present paper.
Another difference between walker trajectories and quantum trajectories is that the quantum description is intrinsically probabilistic and non-classical, while there exist regimes in which the trajectory of the walkers is deterministic and classical (for example, when they bounce exactly in phase with the bath, they can be shown to follow straight lines at constant velocity [14][15][16][17]). However, there also exist regimes in which a Brownian motion is superimposed on their flow lines (e.g., above the Faraday threshold), and other regimes where the trajectories appear to be chaotic [5]. In fact, in several regimes, droplets appear to exhibit ergodic behavior. In practice, ergodicity has been established on the basis of the following observations: if we prepare a walker at the surface of the liquid bath (a corral, for instance), it will progressively explore each part of the surface, following an apparently random motion [4]. If one then visualizes the statistics of the sojourn time of the walker in each of these regions, a striking pattern emerges, bearing more than a simple resemblance to an interference pattern [4,7]. It is this, again remarkable, manifestation of wave-particle duality that first attracted our attention and which lies at the origin of this paper. The onset of quantum equilibrium in the framework of dBB dynamics and in stochastic versions thereof is an important foundational issue in itself, which has motivated numerous studies (see, e.g., [13,[18][19][20][21][22][23][24] as well as [25] and references therein). Several authors in the past have indeed tried to explain how the Born rule emerges from individual trajectories, which is a highly non-trivial problem. In the case of dBB dynamics, it is easy to show that in simple situations the relaxation to the Born statistical distribution does not occur at all, but recent studies [26][27][28][29][30][31] show that in sufficiently complex situations (several modes of different energies for instance) the system might exhibit mixing, which explains the onset of quantum equilibrium in such cases. As we shall show in the present paper, in the case of Nelson-type dynamics, the quantum Brownian motion imposed in such a model accelerates the relaxation to Born's distribution, and in fact ensures that relaxation to the Born rule will almost always occur (as we shall also show). In our view, for the above reasons, de Broglie-Bohm and Nelson-type dynamics are good candidates for explaining how wavelike statistics emerge after averaging a set of apparently chaotic and/or stochastic trajectories.
Briefly summarized, our main goal is to explain the emergence of aforementioned interference patterns in the framework of the dynamical models of de Broglie-Bohm and of a stochastic version thereof which is based on the models of Bohm-Vigier [18] and Bohm-Hiley [19] but which is formally close to Nelson [13]. Both models are introduced in Section 2. Here, it is worth noting that thus far there is no experimental evidence that droplets indeed follow de Broglie-Bohm and/or Nelson trajectories. Our approach therefore differs radically from previous studies on droplets, in the sense that we impose a quantum dynamics by brute force, whereas, until now, the attempt to illustrate how chaos may underlie quantum stochasticity has been a pillar of the research on walkers/droplets. In fact, Nelson's original goal, in proposing his dynamics, was to derive an effective wave equation from the properties of an underlying Brownian motion, as in classical statistical mechanics where a diffusion equation is derived from microscopic properties of the atoms. There actually exists an impressive number of attempts in that direction, as, e.g., stochastic electro-dynamics [5,32,33]. However, there exists (as far as we know) no way to derive an effective Schrödinger equation from hydrodynamical models of droplets.
By choosing exactly the opposite approach, i.e., by imposing quantum-like dynamics on the droplets, we pursue three goals. The first one is to describe the onset of quantum equilibrium (and ergodicity). A second objective is to formulate precise quantitative predictions regarding this relaxation process, which can possibly be validated by future experiments. A third objective is to show, for the first time, that certain dBB trajectories present a deep structural resemblance with certain trajectories that have been reported in the literature for droplets trapped in a harmonic potential.
A short discussion of the onset of equilibrium in de Broglie-Bohm dynamics and the importance of coarse-graining is given in Section 3. In the case of our stochastic, Nelson-type dynamics, we derive in Section 4 a new H-theorem showing the relaxation to quantum equilibrium, which does not rely on coarse-graining and is valid at all scales. We pay particular attention to the ergodicity of trajectories in the case of our stochastic dynamics (which mix properties of the de Broglie-Bohm dynamics with Brownian motion). We apply these ideas to discuss ergodicity in the case of the stochastic treatment of a particle trapped in a harmonic potential (Section 5) and to describe the dynamics of a droplet trapped in a harmonic potential (Section 6). In this latter section (in Section 6.1), we also propose some simple experiments by which one can test the applicability of a Nelson-type dynamics to the context of bouncing oil droplets, and we briefly discuss the problems caused by the presence of zeros in the interference pattern that is encoded in the statistics of the trajectories. In Section 7, we study a situation during which the attractor of the probability distribution is no longer a static eigenstate of the (static) Hamiltonian, and we compare the onset of equilibrium in the dBB and stochastic formalisms in that special framework. In Section 8, we tackle the dynamics of droplets in a 2D harmonic potential through a simple model where the pilot wave is treated as a dynamical object. This constitutes a preliminary attempt, ultimately aimed at establishing a dynamics that would combine stochastic and/or dBB dynamics with a feedback of the trajectory on the wave, which is a fundamental feature of droplet phenomenology that has never been addressed in the framework of dBB or Nelson dynamics. The last section is devoted to conclusions and open questions. A short overview of the numerical methods used in the paper is given in the Appendix A.

The dBB Theory
In the following quick overview of the dBB theory we shall limit ourselves to the case of a single particle. In the dBB theory, particle positions exist at all times and they are merely revealed by position measurements, instead of "originating" with the measurement as the standard interpretation of quantum mechanics would have it. The dynamics is described by a wave function which obeys the Schrödinger equation: where V(x, t) is an external potential and m the mass of the particle, as well as by a position x. In order to reproduce the predictions of standard quantum mechanics, one must have that the positions are distributed according to where P (x, t) is the distribution of particle positions over an ensemble of trajectories. An ensemble satisfying condition (2) is said to be in quantum equilibrium. It is also commonly assumed that (2) is satisfied at some initial time. Therefore, in order to be at (quantum) equilibrium for all t, the condition to enforce is As is well known, the probability density |Ψ(x, t)| 2 satisfies the continuity equation is the (probability) current describing the flow of the probability due to (1). The probability density P, on the other hand, will satisfy a continuity equation where v is the velocity field for the particle. Therefore, Equation The expression (7) for the velocity field is of course not the only possible one: any solution of the form where f is a scalar function, will also give rise to Equation (3) (see [34] for more details). Secondly, if one expresses the wave function in terms of its phase S(x, t) and modulus one finds that (10) and that the velocity of the particle at time t is given by Integrating the system given by Equation (11), we recover the dBB trajectory. From the above, it should be clear that the dBB theory is deterministic. Any stochastic element only comes from our lack of knowledge of the initial positions.
In the context of bouncing droplets, we shall view the external wave generated by the droplet as being in one-to-one correspondence with the "pilot wave" Ψ, which guides the position of the dBB particle.

A Simple Realization of de Broglie's Quantum Thermostat-Nelson Dynamics
As mentioned in the introduction, the trajectories of walkers are often characterized by a non-negligible stochastic (Brownian) component which sets them apart from the smooth dBB trajectories. From this point of view, it seems worthwhile to try to model walkers dynamics in terms of stochastic generalisations of dBB dynamics.
de Broglie himself, in fact, considered such generalizations of the deterministic dBB dynamics (which he called the "quantum thermostat hypothesis") to be highly welcome because they might provide a physically sound picture of the hidden dynamics of static quantum states. For instance, if we consider the position of an electron prepared in the ground state of a hydrogen atom, the dBB dynamics predicts that its position will remain frozen at the same place throughout time, which is counterintuitive to say the least. Adding a stochastic component to its velocity could, in principle, explain why averaging the position of the electron over time is characterized by an exponentially decreasing probability density function, in agreement with the Born rule (provided, of course, that ergodicity is present in the problem in exactly the right proportion). A first proposal in this sense was formulated by Bohm and Vigier in 1954 [18] which, later on, was made more precise by Bohm and Hiley [19], but stochastic derivations of Schrödinger's equation by Nelson [13] (and others [32,33] in the framework of stochastic electrodynamics) can also be considered to provide models of the quantum thermostat. Quoting de Broglie: "...Finally, the particle's motion is the combination of a regular motion defined by the guidance formula, with a random motion of Brownian character... any particle, even isolated, has to be imagined as in continuous "energetic contact" with a hidden medium, which constitutes a concealed thermostat. This hypothesis was brought forward some fifteen years ago by Bohm and Vigier [18], who named this invisible thermostat the "subquantum medium"... If a hidden sub-quantum medium is assumed, knowledge of its nature would seem desirable..." (In [35] Ch.XI: On the necessary introduction of a random element in the double solution theory. The hidden thermostat and the Brownian motion of the particle in its wave.) In this paper, we shall consider a particular model of the quantum thermostat in which, as in the Bohm-Vigier model, a single spinless particle suspended in a Madelung fluid moves with the local velocity of the resulting field, given by Equation (11), and is subjected to fluctuations coming from the latter (cf. Figure 1). However, following Nelson, we shall model these fluctuations by means of a particular stochastic process. To be precise: our model is formally the same as Nelson's in that it relies on the same stochastic process. However, in spirit, it is closer to the Bohm-Hiley model [19] in that we do not assume to be at quantum equilibrium (an assumption which is fundamental to Nelson's theory, as was already pointed out by Bohm and Hiley [19]; see also [20] for a detailed presentation and a comparison of both approaches). This stochastic process is defined on a probabilistic space Ω, characterized by a probability distribution P(x, t) and obeying an Ito stochastic differential equation of the general form: where α is the (constant) diffusion coefficient that characterizes the strength of the random part, and dW (t) is a Wiener process in three dimensions. The function γ(x, t) in (12) is a systematic drift, the so-called osmotic velocity, which we shall fix in the following way. The conservation equation of the probability distribution (which we denote by P, in order to stress the difference with the probability in the dBB dynamics which is denoted by P), can be written as the Fokker-Planck equation: If we now require that the quantum equilibrium P(q, t) = |Ψ(q, t)| 2 be a solution of this Fokker-Planck equation, we obtain from Equations (4), (10), and (13) that which is a constraint on the osmotic velocity. The simplest solution of this constraint is [36] γ(x, t) = α 2 In the rest of the paper, we choose the osmotic drift velocity to be Equation (15), with α an a priori free parameter, Nelson's choice for α (α =h/m) being irrelevant if we apply this formalism to droplets.
In summary, our Nelson dynamics is fully defined by the following Ito equation where dW i (t) represents a Wiener process with and by the Fokker-Planck equation for the associated probability distribution P(x, t) where Ψ(x, t) satisfies the Schrodinger equation: At quantum equilibrium, i.e., when P(x, t) = |Ψ(x, t)| 2 , the diffusion velocity is balanced by the osmotic term and the Bohm velocity is recovered, on average.
We shall now discuss the details of the relaxation towards quantum equilibrium, in the dBB and stochastic formalisms.

Relaxation to Quantum Equilibrium in the de Broglie-Bohm Theory
In our presentation of the dBB theory for a single particle, in Section 2.1, we assumed that the particle positions are initially distributed according to Born's law over an ensemble. The dynamics then ensure that the same relation will hold for any later time. This is the assumption de Broglie and Bohm made in their original papers [11,12,37]. Although Bohm tried, already in the 1950s (first on his own-see, e.g., [11] (Section 9)-and later with Vigier [18]), to relax this assumption by modifying the dynamics, to many authors working today on the dBB theory it is still an assumption which has to be made (the final objective of de Broglie, Bohm, Vigier, and Nelson-and many other contributors to various realistic hidden variable interpretations in which quantum systems are assumed to be localized in space at any time-was to rationalize wave-like statistics in terms of individual trajectories; the same problem occurs in droplets phenomenology and, according to us, admits no fully satisfying solution yet). According to Valentini [21,38,39], however, there is no need to assume that the particle positions are initially distributed according to Born's law or to modify the dynamics. His claim is that an ensemble in which Born's law is not satisfied (so-called quantum non-equilibrium) will evolve naturally towards quantum equilibrium, provided that the wave function leads to sufficiently complex dynamics. This relaxation process has to take place on a coarse-grained level and can only occur if the initial distributions do not display any fine-grained micro structure.
Let us first explain the need for coarse-graining. Let us introduce the function f = P /|Ψ| 2 , as in [26]. An important implication of (6) is that the function f is conserved along the dBB trajectories: Hence we have that where x i is the initial position of the particle which leads to x, when evolving from t i to t according to the dBB dynamics. If one assumes that P (x i , t i )/|Ψ(x i , t i )| 2 = 1, relaxation to quantum equilibrium is clearly impossible, at least at the microscopic level. However, as argued by Valentini [21], relaxation is possible at the coarse-grained level, provided the initial distribution does not display any fine-grained microstructure. The operational definition of the coarse-graining is as follows. We divide the domain of interest A ⊂ Ω into small cubes of equal edge length (we call them coarse-graining cells, or CG cells for short). These CG cells do not overlap and their union is equal to A. The coarse-grained densities, which we denote by P (x, t) and |Ψ(x, t)| 2 , are then defined as where the domain of integration is the CG cell containing x.
We can now discuss the second assumption: the lack of a fine-grained microstructure in the initial distribution. Let us assume we have a non-equilibrium distribution P (x, t i ) which relaxes to quantum equilibrium at the coarse-grained level, under the dynamics generated by the wave function Ψ(x, t). As the dBB theory is time-reversal invariant, in the time-reversed situation, under the dynamics generated by Ψ * (x, −t), we would have a distribution that moves away from quantum equilibrium. Thus, it would seem that time-reversal invariance contradicts the possibility of relaxation to quantum equilibrium. This conclusion is unwarranted, however: as the initial distribution P (x, t i ) relaxes to quantum equilibrium, it retains information on the original values of f (which are constant in time) and thereby acquires a fine-grained microstructure, which means that at the final time t f , P (x, t f ) will differ significantly from P (x, t f ). Therefore, in the time-reversed situation, the initial distribution would exhibit a fine-grained microstructure, which is prohibited under our assumption, thereby breaking the time-reversal invariance.
In order to quantify the difference between the distribution P (x, t) and the quantum equilibrium condition |Ψ(x, t)| 2 at the coarse-grained level, Valentini [21,38,39] introduced the entropy-like function where P and |Ψ| 2 as in Equations (23) and (24), for which he has shown the (quantum) H-theorem, under the assumption of no fine-grained microstructure. It should be stressed, however, that this is not necessarily a monotonic decay and therefore does not prove that quantum equilibrium will always be reached. It merely indicates a tendency towards relaxation. The strongest support for the idea of relaxation to quantum equilibrium comes from numerical simulations of the evolution of non-equilibrium distributions for various quantum systems [26][27][28][29][30][31] (see [25] and references therein for a review). The first numerical simulations were performed by Valentini and Westman [26] who showed, in the case of a 2D box, that relaxation quickly takes place for a wave function which is a superposition of the first 16 modes of energy (the superposition being equally weighted). It was also hinted that the nodes of the wave function, with their associated vorticity, play a crucial role in the relaxation process, as purveyors of chaos (or mixing) in the dynamics. This later claim was properly understood in [40]. The dependence of the relaxation timescale on the coarse-graining length and on the number of energy modes was studied in [28]. In [31], it was shown that quantum systems with a low number of modes are likely to never fully relax, in which case H reaches a non-zero residue value. However, such a scenario becomes unlikely as the number of modes increases. According to the quantum non-equilibrium hypothesis, standard quantum mechanics is only one facet of the pilot-wave theory, that of quantum equilibrium, leaving the possibility for new physics: that of quantum non-equilibrium. One should assume of course that during our time we have only had (or can only have) access to systems for which quantum equilibrium has already been reached. But that does not mean that quantum non-equilibrium never existed in the early universe (which could be inferred from the observation of the remnants of the early fractions of seconds of the universe, just after the Big Bang [41]), or that some, yet undetetected, exotic quantum systems cannot still be in quantum non-equilibrium today [42]. This is why droplets are appealing, because if their dynamics do present analogies with dBB dynamics, their study will allow us to observe relaxation to (quantum) equilibrium "in real time" in the lab with the naked eye, which is not possible with quantum systems for which we have no direct access to individual trajectories.

An H-Theorem for Nelson Dynamics
Let us start by introducing an analog of Valentini's entropy, Equation (25), for the probability distribution P(x, t) associated with our Nelson dynamics, as defined by Equations (16)(17)(18)(19): which is a special instance of a relative entropy known as the Kullback-Leibler divergence [43]. We also define a second non-negative functional where Note that we always impose the boundary conditions |Ψ| 2 ∂Ω = P ∂Ω = 0 and f ∂Ω = 1 so as to avoid divergence of these integrals on the boundary of Ω.
It should be noted that the entropy of Equation (27) or the functional of Equation (28) we shall use to quantify the relaxation to quantum equilibrium are very different from the entropies usually considered in the context of classical H-theorems (like e.g., the Boltzmann entropy). One should bear in mind, however, that quantum equilibrium is radically different from classical equilibrium [44] and has no connection whatsoever with relaxation to quantum thermal equilibrium, for the simple reason that the Born distribution of positions reached by an ensemble of trajectories à la Nelson or dBB is not a thermal distribution.
To understand why the functionals in Equations (27) and (28) are non-negative and why they are zero if and only if (quantum) equilibrium is reached (that is to say when f = 1 everywhere in space), it is important to note that the integrands of H V and L f satisfy the inequalities for which any of the possible equalities only hold when P = |Ψ| 2 . (This is immediate from the trivial inequality: Now, since both P(x, t) and |Ψ(x, t)| 2 are probability distributions, i.e., since we have Ω Pdx = Ω |Ψ| 2 dx = 1, it follows from (30) that whenever H V (t) and L f (t) are well-defined, they satisfy the following inequalities: Moreover, for the same reason, L f can be re-expressed as Ω d 3 x f (P − |Ψ| 2 ) − (P − |Ψ| 2 ) , the integrand in which is non-negative due to Equation (30). Therefore, L f can only be zero if its integrand is zero, i.e., if P = |Ψ| 2 (if P, |Ψ| 2 , and f are sufficiently smooth, which is something we shall always assume unless otherwise stated). Similarly [21], one also has that H V can only be zero when P = |Ψ| 2 everywhere in Ω.
Let us now prove the relaxation to quantum equilibrium. Substituting P = f |Ψ| 2 in the Fokker-Planck Equation (18), and using Equations (4) and (10), it is easily verified that Rewriting L f as its behavior in time can be calculated using Equations (32), (4), and (10): which is of course strictly negative, for all t, as long as ∇ f and |Ψ| 2 are not identically zero. Hence, if |Ψ| 2 is not zero throughout Ω, L f will decrease monotonically for as long as f is not (identically) equal to 1 on Ω, and therefore necessarily converges to 0, a value it can only attain when f ≡ 1 or, equivalently, when P ≡ |Ψ| 2 . We have thus established a strong H-theorem showing that, in the case of Nelson dynamics, any probability distribution P necessarily converges to |Ψ| 2 , if the latter does not become zero identically. Note that this excludes the case of a free particle for which lim t→+∞ |Ψ(x, t)| 2 = 0, for all x, which means that dL f dt tends to zero even when f does not converge to 1. A result, similar to the above, is also easily established for H V since L f dominates the latter, or alternatively from the formula The above results show that (excluding the case of the free particle) Nelson dynamics, naturally, exhibits relaxation towards quantum equilibrium and that it does so for general initial probability distributions (at least, as long as the initial distribution is smooth enough). In this stochastic setting, there is therefore no need for any assumptions on the microstructure of the initial distributions, nor is there any need for the coarse-grained hypothesis when deriving an H-theorem.
Note that these results also show that we have, in fact, convergence of the distribution P to the quantum equilibrium distribution |Ψ| 2 in the L 1 norm. This is a consequence of the so-called Csiszár-Kullback-Pinsker inequality [43]: where This generalizes the results by Petroni and Guerra [22,23] obtained in their study of the relaxation towards quantum equilibrium in the framework of the Nelson dynamics of a single particle in a harmonic potential. The L 1 norm is also used by Efthymiopoulos et al. [25] in the context of the dBB theory.
We shall illustrate these results by means of numerical simulations for the case of a ground state for the 1D-harmonic oscillator in Section 5.4, for the case of the 2D-harmonic oscillator in Section 6.2, and in the case of a coherent state in Section 7.1.
A last important remark concerns the influence of possible zeros in the equilibrium distribution Ψ(x, t), which would give rise to singularities in the osmotic velocity terms in the Ito equation, Equation (16), or in the Fokker-Planck equation, Equation (18) (or equivalently in Equation (32)), and might make the functions H V and L f ill-defined. In Section 6.2, we discuss the case of the first excited state of the 1D-harmonic oscillator, for which Ψ(x, t) has a node at x = 0. One could in fact imagine studying higher excited states for which one has a finite number of nodes. In that case, the osmotic velocity of Equation (15) will have simple poles at a finite number of positions in x. At the level of the Ito equation, one would not expect a finite set of poles to cause any particular problems, not only because the probability of hitting a pole exactly in the stochastic evolution is zero but also because the osmotic term tends to move the particle away from the pole very quickly. Similarly, a finite number of simple poles in the convection-diffusion equation, Equation (32), for f only influence the velocity field in the convection term in a finite number of distinct places, and it is to be expected that this would have the effect of actually enhancing the mixing of information in the system.
Moreover, it is also clear that simple nodes in Ψ(x, t) only give rise to (a finite number of) logarithmic singularities in the integrand of H V and that the integral in Equation (27) therefore still converges. The H-theorem for H V derived above is thus still valid, and an arbitrary distribution P (sufficiently smooth) will still converge to quantum equilibrium, even in the presence of nodes for Ψ(x, t). The same cannot be said, however, of the function L f , as simple zeros in Ψ(x, t) give rise to double poles in the integrand and a possible divergence of the integral of Equation (28). Hence, at the beginning of the evolution, for an arbitrary P, the function L f might take an infinitely large value (the integrand only diverges when |Ψ| 2 P, i.e., when it is positive) but as soon as convergence sets in (which is guaranteed by the H-theorem for H V ), the divergent parts in its integrand will be smoothed out and the function L f will take finite values that converge to zero as time goes on. Of course, when calculating these quantities for the results of numerical simulations, there is always some amount of coarse-graining going on and genuine infinities never occur.

Relaxation to Quantum Equilibrium and Nelson Dynamics: The Static Case
In this section, in order to simplify the discussion, we will only consider the case of stationary states Ψ st (x) for the one dimensional Schrödinger equation, i.e., energy levels for which S = −E t and which therefore have zero Bohm velocity (11): ∇S ≡ S x = 0.

Fokker-Planck Operator and a Formal Connection to a Schrödinger Equation
There exists a wide literature [45,46] concerning a particular method for studying the convergence of solutions of the Fokker-Planck equation to a stationary one, which is only sporadically mentioned in the literature devoted to Nelson dynamics [47]. This approach makes it possible to quantify very precisely the speed of convergence to equilibrium, in terms of (negative) eigenvalues of the Fokker-Planck operator. In order to show this, let us rewrite the Fokker-Planck Equation (18) in terms of the Fokker-Planck operator L: where (15): Note that, due to the presence of the first derivative ∂ ∂x , the L operator is not Hermitian. Now, in order to establish the H-theorem, we must prove that in the long-time limit this equation tends to a stationary solution P st = |Ψ st | 2 . The key idea here is to transform the Fokker-Planck equation to a simple diffusion equation through the transformation under which the r.h.s. of Equation (40) reduces to where H st is now a Hermitian operator: The function g(x, t) thus obeys a "Schrödinger-like" equation (though with imaginary time) with an effective potential ( H st ) that depends on γ(x): Note that the effective potential is exactly the Bohm-quantum potential defined by which can be expressed in terms of the osmotic velocity (41) as

Superposition Ansatz
We can now represent the solution of Equation (45) as a superposition of discrete eigenvectors (all orthogonal, as the operator H st is Hermitian) and impose the superposition ansatz [48]: Equation (45) is separable and gives rise to the eigenvalue problem: As a result, we have for a set of constants a k and where all the λ k are real (as H is Hermitian), for eigenfunctions g k (x) that satisfy the orthonormality conditions: Thus, we have the expression By construction, the function P st (x) is an eigenstate of the effective Hamiltonian with energy 0. We shall associate the label λ 0 with this energy level.
In order to have a well defined probability distribution and to avoid any divergence in time, it is clear that all eigenvalues −λ k have to be negative, which requires Ψ st to be the ground state of the effective Hamiltonian H st : just as in the case of the usual Schrödinger equation, the eigenvalues −λ k in Equation (49) are all negative only if Ψ st (x) has no zeros (see also [49], Appendix 2, for an elementary proof that all λ k are indeed positive if Ψ st (x) h as no zeros).
If Ψ st (x) does have zeros, the osmotic velocity will have singularities. In [49] (Appendix 2), we consider what happens in the case when Ψ st (x) is an excited state of the harmonic oscillator and we derive a formal solution in terms of the eigenvalues −λ k , which are now not all negative, thus revealing the appearance of instabilities for cases where the above formalism would still be valid.

One-Dimensional Oscillator and the Evolution of Gaussian Distributions for the Ground State
In [49] (Appendix 2), we discuss the application of the method of the effective Hamiltonian outlined in Section 5.1 to this particular problem, and we derive a Green function for the associated Fokker-Planck equation when Ψ st is the ground state of the one-dimensional oscillator. This Green function K P (x, x , t) is defined through where the kernel K P is given by An important property of the Green function for this case is that, if |Ψ(x)| 2 and P(x, 0) are Gaussian, then P(x, t) will still be Gaussian (53). Let us define the ground state as for which we can then write Injecting Equation (56) in the Fokker-Planck equation, Equation (40), gives a differential equation for which is readily solved as as well as an equation for b(t), From Equations (56) and (60), we can then calculate the width of the non-equilibrium Gaussian as where σ 2 eq represents the width 1/(4a) of the equilibrium distribution of Equation (55). Clearly, x t→∞ = x eq = 0 with a characteristic relaxation time inversely proportional to the diffusion coefficient α. Moreover, which has the same sign as that of the difference (σ eq − σ x (0)). Hence, σ x (t) converges monotonically to the equilibrium value σ eq , with a characteristic time inversely proportional to the diffusion coefficient α, as can be seen in Figure 2.  (16), for the ground state (Equation (55)) of the 1D harmonic oscillator), whose initial positions are normally distributed, for 5 different choices of distribution width (for a = 0.5 and α = 1). We observe, in each case, convergence to the equilibrium of Equation (55) as predicted by the theory.

Ergodicity in the Relaxation to Quantum Equilibrium for the Ground State of the Harmonic Oscillator
We have just shown how Gaussian initial distributions converge towards quantum equilibrium, but one could also ask the same question for non-Gaussian initial distributions. Convergence is guaranteed by the H-theorem of Section 4, but contrary to the Gaussian case, we have no clear measure for the rate of convergence, except for the entropy-like functions H V (27) and L f (28), or the L 1 norm (39), defined in Section 4. The evolution, in time, of these three quantities is shown in Figure 3, for the stochastic trajectories obtained from 20,000 uniformly distributed initial conditions. The relaxation towards quantum equilibrium is clearly visible in all three quantities. As expected, the convergence of H V is extremely fast. Note that, although initially very large, L f quickly matches L 1 , up to numerical fluctuations.  One important question concerning this relaxation process is of course that of possible ergodicity. Since we want to study the ergodic properties of Nelson dynamics in a numerical way, we choose the definition of ergodicity that is, in our approach, the easiest to test. Defining the time average h of a function h on Ω by the limit (if it exists), (where x t represents the position of a particle at time t , as obtained form the Ito stochastic differential equation, Equation (16), for an initial condition x) we say [50] that the corresponding stochastic process is ergodic if the time average of any bounded function h on Ω is always independent of x.
Since for bounded h the time average is also invariant under shifts in time, we can say that we have ergodicity if all time averages of such functions are in fact constants. The main reason for choosing this particular definition is that it is well-suited to empirical testing, since it is of course sufficient to establish constancy of the time averages for all indicator functions χ A of arbitrary (measurable) sets A ⊂ Ω, for the analogous property to ensue automatically for all bounded functions on Ω. Another reason for choosing this particular definition is that it is also applicable to non-stationary stochastic processes, as in the case of the coherent state of Section 7.
More precisely, we need to verify that is independent of both t and x, for any measurable A ⊂ Ω. Remember that one has of course that In the present case, i.e., that of the Nelson dynamics defined by the stationary (ground) state of the 1D harmonic oscillator, it is clear that the distribution |Ψ st | 2 obtained from the ground state eigenfunction Ψ st is a stationary solution to the associated Fokker-Planck equation, Equation (18). This distribution provides a natural invariant measure µ on Ω: dµ = |Ψ st | 2 dx, for which Ω dµ = 1 and If a stationary stochastic process is ergodic, i.e., if for such a process all χ A are indeed constants, the values of these constants are simply the measures of the subsets A [51]. Therefore, when one needs to decide whether or not a stationary stochastic process is ergodic, it suffices to verify that χ A = µ(A), for any A ∈ Ω.
The usual way to check this condition is to consider sampling time averages for a sufficiently refined "binning" of Ω. Starting from a particular initial particle position x, we calculate the trajectory x t that follows from the Ito stochastic equation, Equation (16), for a sufficiently long time t. As was explained for the coarse-graining in Section 3, the configuration space Ω is subdivided into a large number of non-overlapping cells or "bins" A k (k = 1, . . . , N b ), each with the same volume ∆x.
The trajectory x t (t ∈ [0, t]) is then sampled at regular intervals ∆t, yielding N + 1 sample positions x n∆t (n = 0, . . . , N), for N = t/∆t. We then define the sampling function ϕ N,k which is a discretization of 1 t t 0 χ A x t dt in Equation (64) and which gives the frequency with which the (sample of the) orbit visited the bin A k . Hence, if in the limit N → +∞, for diminishing bin sizes ∆x and sampling steps ∆t, the normalized distribution obtained from ϕ N,k /∆x tends to a constant distribution (and, in particular, does not depend on the initial positions x), then the stochastic process is ergodic according to the above definition.
Moreover, since in that case χ A k = µ(A k ), this normalized distribution must in fact coincide with that for the invariant measure for the stationary process. For example, in the case at hand, if the normalized distribution we obtain is indeed independent of the initial positions, then since µ(A k ) = |Ψ st (x)| 2 x=ξ ∆x for some point ξ ∈ A k , we must have that for sufficiently large N i.e., the empirical distribution obtained from this sampling time average must coincide with the stationary quantum probability |Ψ st | 2 . This is exactly what we obtain from our numerical simulations, as can be seen from the histograms depicted in Figure 4. After a certain amount of time, the histograms we obtain indeed converge to the equilibrium distribution, and this for arbitrary initial positions. The convergence clearly improves if we increase the integration time, or if we diminish the spatial size of the bins (while diminishing the sampling time step in order to keep the occupancy rate of each bin high enough). Although purely numerical, we believe this offers conclusive proof for the ergodicity of the Nelson dynamics associated with the ground state of the harmonic oscillator in one dimension.  The same can be said, in fact, for the 2-dimensional oscillator which will be the main topic of Section 6 (and even for the 2D corral as can be seen in [49]). Some results of a simulation of a single trajectory under the Nelson dynamics for the ground state of this system are shown in Figure 5, in which the red dot in the plot on the left-hand side indicates the (final) position of the particle at time t. The probability distribution obtained by sampling the trajectory clearly decreases with the distance to the origin along concentric circles.

2D Harmonic Oscillator
Experimentally, it has proven possible to study the dynamics of bouncing droplets under the influence of an effective harmonic potential in two dimensions, thanks to a well-chosen electro-magnetic configuration and magnetic droplets [52]. It is therefore interesting to compare predictions that we, on our side, can make in the framework of Nelson dynamics, with actual experimental observations of droplets dynamics (see [53] for a pioneering work very similar to ours in the case of the double slit experiment). We think that an important comparison to make concerns the convergence to equilibrium.
For example, if the initial distribution of positions projected along a reference axis, say X, fits a mixture of the ground state and the nth Fock state (n = 1, 2 · · · ) (Appendix 2, [49]) for the 2D harmonic oscillator (conveniently weighted in order to respect the ineluctable constraint of positivity), our Nelson-like model predicts that the typical time of convergence to equilibrium will scale like the inverse of the eigenvalue of the nth Fock state, i.e., as 1/n, which constitutes a very precise quantitative prediction. This follows from Equation (52), when P st (x) is the Gaussian ground state of the 1D harmonic oscillator and where the eigenfunctions g k are the energy eigenstates (Fock states) of the harmonic oscillator (this, of course, because of the separability of the Schrödinger equation and of our Nelson dynamics along X and Y in the case of an isotropic 2D oscillator).
A possible way to measure this characteristic time would be to record the projections along X of trajectories that correspond to an equally spaced grid of initial positions, weighted so as to fit a mixture of the ground state with the nth Fock state (n = 1, 2 · · · ), and to compare the histogram constructed in this way at different times with theoretical predictions derived from Equation (52).
Another precise quantitative (theoretical) prediction, which is even simpler to verify, is that if we prepare a droplet many times at exactly the same initial position, the position obtained after averaging over all trajectories will (1) decrease exponentially in time and (2) be characterized by a decay time which scales like 1/aα, by virtue of the discussion in Section 5.3 and particularly Equation (58). It has been suggested that droplet trajectories might be characterized by a quantum-like Zitterbewegung, which can be seen in relativistic quantum dynamics as an intrinsic Brownian motion at the Compton scale [54,55], and various proposals have been formulated in order to express the amplitude and frequency of this Zitterbewegung [5,56] in terms of the parameters characterizing droplet dynamics (these are, e.g., the viscosity of the fluid, the mass of the droplets, the ratio between the amplitude of the vibrations imposed on the bath and the Faraday threshold, the oil temperature [5,57], and so on). Exploring these analogies in depth lies beyond the scope of this paper, but the aforementioned attempts ( [56] in particular) pave the way for introducing a Brownian component in the description of droplet trajectories.

Presence of Zeros in the Interference Pattern
One of our first motivations, when we decided to incorporate a Brownian component in the dBB theory in order to simulate the dynamics of droplets, was the pioneering paper [4] reporting on observations of a walker trapped in a spherical 2D cavity (corral), for which the histogram of positions occupied over time by a single droplet trajectory faithfully reproduces the Bessel function J 0 (this is also related to the Green function of the Helmholtz equation, with a typical length equal to the Faraday wave length of the vibrating bath over which droplets propagate [16]). These observations reveal, in a telling way, the presence of a pilot-wave that guides the dynamics of the particles, and raise the question of ergodicity.
If we try the approach we used for the 2D harmonic oscillator in the case of the corral (effectively replacing the Gaussian ground state of the 2D harmonic oscillator by the zero order Bessel function), we are immediately confronted with problems caused by the presence of zeros in the Bessel function. In particular, the eigenvalues −λ k of the Fokker-Planck operator in Equation (49) are not always negative when zeros are present, which of course would menace the stability of the relaxation process. However, as we already indicated in Section 4, although the effect of zeros of the pilot wave in our Nelson dynamics is by no means trivial, there are several observations that indicate that this problem is not crucial.
First of all, as mentioned in Section 4, the Wiener process makes it in principle possible to "jump" over the zeros of the equilibrium distribution. This has actually been confirmed in numerical simulations for the case of the 1D harmonic oscillator, where we imposed that the equilibrium distribution P st is the square modulus of the first excited (Fock) state ( [49], Appendix 2A), with the following amplitude: Indeed, as can be clearly seen from Figure 6, the particle will, from time to time, jump over the zero in the middle (with the same probability from left to right as in the opposite direction), in such a way that finally the trajectory covers the full real axis, while the histogram of positions faithfully reproduces the quantum prediction P st = |Ψ st | 2 = |Ψ 1 (x, t)| 2 . This indicates that, even in the presence of a zero in the equilibrium distribution, the relaxation process is still ergodic. The relaxation of a uniform initial distribution to this quantum equilibrium is shown in Figure 7, for the quantities H V , L f , and L 1 . A second indication that the problem posed by the presence of zeros is not so serious stems in fact from the experimental observations. Indeed, if we study the observations reported in [4] for the case of a corral, it is clear that the minima of the histogram expressing the distribution of positions of the droplet are in fact not zeros. This, undoubtedly, is due to the presence of a non-negligible residual background. Without this background, the droplet would never pass between regions separated by zeros: due to the rotational symmetry of the corral, the zeros form circles centered at the origin and the position histogram obtained from a trajectory would remain confined to a torus comprising the initial position. This, however, is clearly not the case, which thus suggests the following strategy: to simulate Nelson dynamics with a static distribution P st = |Ψ st | 2 given by the Bessel function J 0 but supplemented with a constant positive background : In this case, the singularities of the Fokker-Planck equation automatically disappear and, despite the fact that we have no analytic expression for the solutions as in the case of the ground state of the harmonic oscillator, we are able to numerically simulate Nelson dynamics without difficulty. The results of these simulations are shown in Figure 8. The osmotic velocity in the Nelson dynamics clearly tends to bring back the particle to regions where |Ψ| 2 has extrema and the resemblance with the plot on the left is striking. The fact that this result again does not depend on the choice of initial condition strongly suggests that the relaxation process to quantum equilibrium is also ergodic in this case.

Nelson Dynamics and Asymptotic Coherent States
Up to now, we have developed analytic and numerical tools aimed at studying the onset of equilibrium when the asymptotic equilibrium distribution is static. As the H-theorem of Section 4 is also valid for non-stationary processes, one of course expects relaxation to take place even if the asymptotic state is not static, for instance, when it is a Gaussian distribution, the center of which periodically oscillates at the classical frequency ω of the oscillator without deformation (typical for coherent states). In fact, our numerical simulations show not only that equilibrium is reached even in this case, but also that this relaxation is ergodic.
More precisely, we considered a wave function in the coherent state where ϕ is a global phase containing the energy, andx t (p t ) is the mean position (momentum) of a classical oscillator at time t:x t =x 0 cos (ωt) andp t = −mx 0 sin (ωt) , with ω = 2aα (α =h/m). For this ansatz we solved the Ito equation, Equation (16), numerically for a collection of initial conditions. As can be seen in Figure 9, the trajectories are affected by the stochastic evolution but keep oscillating at the same period because of the deterministic part of the Ito process. Notice, however, that the trajectories seem to be approaching classical trajectories that only differ from each other by a simple shift. This can be explained as follows: at equilibrium (cf. Figure 10), the Brownian motion is balanced by the osmotic velocity and the dBB velocity is recovered "on average." Now, the center of the Gaussian distribution moves at a classical velocity by virtue of Ehrenfest's theorem; moreover, in the present case, the dBB velocities can only depend on time and not on space as the envelope of a coherent state moves without deformation. Hence, the dBB trajectories obtained at equilibrium are, in fact, classical trajectories that only differ by a mere shift in space (the magnitude of which, however, may change over time). x(t) Figure 9. Numerical solutions of the Ito stochastic differential equation, Equation (16), corresponding to the coherent state of Equation (70), for three different initial conditions. We usedx 0 = 1, a = 0.5, and α = 1 and expressed the results in natural units.
Secondly, as can be clearly seen on Figure 10, even for a uniform initial probability distribution, the convergence to the quantum equilibrium is remarkably fast and the converged distribution faithfully follows the oscillating motion of the non-stationary equilibrium distribution. The remarkable speed of the convergence to quantum equilibrium is corroborated by the decay of the functions H V and L f and of the L 1 norm shown in Figure 11. Moreover, Figure 12 depicts the sampling time average (as defined in Section 5.4) of a single trajectory for this non-stationary stochastic process. The convergence of the sampling distribution to a static distribution Φ(x), described by the integral of |Ψ(x, t)| 2 as given by Equation (70), over a period of the oscillation is striking. As the asymptotic distribution Φ(x) does not depend on the choice of initial condition, we conclude that the relaxation to equilibrium for the non-stationary stochastic process associated with Nelson dynamics for the coherent state (70) is ergodic (in the sense explained in Section 5.4).

Onset of Equilibrium with a Dynamical Attractor in dBB Dynamics and Nelson Dynamics
If one wants to investigate the onset of equilibrium in dBB dynamics, one obviously has to consider non-static asymptotic distributions since in static cases the dBB dynamics freezes the trajectories (as the phase of the wave function is then position-independent). Even in the case of a coherent state (see Section 7.1), the distribution of dBB positions would never reach equilibrium because it moves as a whole (as the shape of a coherent state remains the same throughout time). In a sense, coherent states behave as solitary waves. Moreover, the absence of zeros in the wave function might explain why mixing does not occur. In Figure 13 we show the result of simulations of dBB trajectories in the case of a 2D harmonic oscillator for a quantum state consisting of a superposition of equally weighted products of states along X and Y, chosen among M energy (Fock) states ( [49], Appendix 2A), with randomly chosen initial phases θ n x ,n y : e i θ nx ,ny −iω(n x +n y +1) t ψ n x (x) ψ n y (y) . We then compared the relaxation process for dBB with the quantum thermostat given by Nelson dynamics for M = 4 2 = 16 energy states. The results are shown in Figure 14 in which H V (for the dBB and for the Nelson dynamics) and L 1 (for both the dBB and Nelson dynamics) are plotted at the (same) coarse-grained level. We started from a uniform distribution of positions; we took α = 0.1. In both cases, the position distributions P and P converge to |Ψ| 2 . Moreover, we recover an exponential decay for H V , as already observed in [26], even in the absence of stochastic (Brownian) noise à la Nelson. However, we observe that the convergence to equilibrium occurs faster in the presence of the quantum thermostat.

Dynamical Model for Droplets and Double Quantization of the 2-D Harmonic Oscillator
In this section we shall focus on the description of droplets dynamics as described in [9,52], for a magnetized droplet moving in an isotropic 2-D harmonic potential. We shall show that dBB dynamics allows us to reproduce some of the main features of the experimental observations. In [9,52], it is reported that stable structures appear in the droplets dynamics whenever a double quantisation condition is satisfied. The Hamiltonian of the isotropic 2-D harmonic oscillator being invariant under rotations, we may indeed impose a double quantisation constraint, requiring that the energy states of the 2D quantum harmonic oscillator are also eigenstates of the angular momentum. In polar coordinates, these states (which are parameterized by two quantum numbers, the energy number n and the magnetic number m) are expressed as follows [58]: where L |m| k are the generalized Laguerre polynomials and k = n−|m| 2 . Note that these solutions are linear combinations of the product of Fock states in x and y.
A first experimental result reported in [9] is the following: trajectories are chaotic and nearly unpredictable unless the spring constant of the harmonic potential takes quantized values that are strongly reminiscent of energy quantization (under the condition that, during the experiment, the size of the orbits is fixed once and for all). For quantized energies-in our case given by E n = (n + 1)hω, for some "effective" value ofh to be determined from actual experiments-stable orbits, to which one can attribute yet another quantum number, appear, this time for the angular momentum, which is strongly reminiscent of the magnetic number (the eigenvalue of the orbital momentum, perpendicular to the surface of the vessel, is given by the product ofh and m). In [9] it is shown, for instance, that for the first excitation (n = 1, m = ±1) droplet orbits are circular or oval, turning clockwise or anti-clockwise depending on the sign of m. At the second energy level (n = 2, m = −2, 0, +2), ovals appear again for m = ±2 and lemniscates for an average value of the angular momentum < m > = 0. At the fourth energy level (n = 4, m = −4, −2, 0, 2, 4), trefoils appear (for m = ±2).
We simulated dBB trajectories, always considering a superposition of one of the aforementioned doubly quantized eigenstates ψ n,m with the ground state: where ϕ j and ξ j are real numbers with 0 < ξ 0 ξ j =0 . Computing the guidance relation of Equation (11) for a single eigenstate (74), one ends up with a value for ∇S for which the trajectories are circles of radius R around the origin, with tangential velocities proportional to m/R. In particular, the dynamics is frozen when m = 0.
Mixing the wave function with the ground state, however, generates a periodic (in time) component in the dBB velocity field, which turns circular orbits into ovals when ξ 0 is small enough, and eventually generates more complex structures, such as "rosaces" instead. We also tuned the energy difference between the ground state and the excited states such that two timescales characterize the dynamics. These are the "centrifugal" period, necessary for drawing a full circle around the origin, which varies as m/R 2 , and the "Bohr" period which varies like T/(n + 1), where T is the classical period of the oscillator. Tuning these parameters, we were able to simulate dBB trajectories very similar to those reported in [9]. For instance, we found circles and ovals (see Figure 15a,b) for (n, m) = (1, 1) and (n, m) = (2, 2). Note that the lemniscate cannot be obtained with a superposition of the ground state and the (n, m) = (2, 0) state for which dBB velocities are necessarily purely radial, contrary to the suggestion made in [9]. Instead, it should be generated with a superposition of the ground state with (n, m) = (2, +2), (2, −2), and (2, 0) in which the weights of the m = +2 and −2 components are slightly different (see Figure 15c). Figure 16 shows further detail of the evolution along this trajectory. Tuning the energy, we were also able to generate a trefoil and a "rosace" (see Figure 17). Each plot is associated with a different combination (n, m), as indicated. In the (a,b) graphs, we imposed a = 1 and, respectively, ω = 1, ξ 0 ξ 2 = 0.05 and ω = 0.5, ξ 0 ξ 3 = 0.05; for (c), we imposed a = 3, ω = 0.5, ξ 0 ξ 3 = 0.0708, ξ 0 ξ 2 = 0.0456, and ξ 0 ξ 1 = 0.0773. It is worth noting, however, that chaos is omnipresent in the dBB dynamics for this system, in the sense that the trajectories exhibit an extreme sensitivity to the initial conditions, which explains why these dBB orbits mimicking stable droplets orbits are in general unstable. For instance, Figure 18 shows intermittent transitions between an oval trajectory and a lemniscate (as has also been reported in [9]), for a superposition of the ground state with the (n, m) = (2, +2), (2, −2), and (2, 0) states. Preliminary results furthermore show that the trajectories are also unstable under Nelson dynamics, i.e., in the presence of "noise," whenever this noise (parameterized by α in (16)) exceeds a critical value. Note that many experiments involving droplets are characterized by a lack of stability and predictability. For instance, the appearance of interferences similar to those obtained in a double slit experiment (see [49,59] for a description à la Nelson of the double slit experiment) has been attributed to "air currents" in [60]. Therefore, although our approach might not explain every detail of the double quantization reported in [9], it does reproduce many of its essential features, and we believe it would be very interesting to deepen this analogy. For instance, having access to the empirical values of the weights of the ground state, or of the effective values ofh and of the mass in the case of droplets [56], would allow us to test our model in real detail. Another experiment, reported in [61], during which both the position of the droplet and the excitation of the bath are monitored, and where a superposition between two distinct modes of the bath is reported, could also provide more insight and might offer some means to test the validity of our model: using exactly the same observation device, but this time in the case where the droplet undergoes a 2-D isotropic potential, would allow one to check whether the modes of the bath are similar to the (n, m) quantum modes which we associate with the quantized droplets trajectories.

Conclusions and Open Questions
In this paper we studied stochastic, Nelson-like dynamics and dBB dynamics, with the aim of simulating the dynamics of droplets. The stochastic approach has the merit that it explicitly takes into account the influence of noise on the dynamics [59,62]. In contrast to experiments where noise is considered to be a parameter that should be minimized, here, noise is considered to be a relevant parameter for the dynamics (see also [53]). For instance, as we have shown, it plays an essential role in the relaxation towards equilibrium and in the ergodicity of the dynamics. In the dBB approach, on the other hand, the main ingredient is the chaotic nature of the dynamics [25]. Both models thus shed a different light on the dynamics and could possibly fit diverse sets of regimes in droplet dynamics. Note that in the limit where the amplitude of the Brownian motion in our Nelson dynamics tends to zero, the dynamics approaches dBB dynamics very closely. In sufficiently complex situations (e.g., when the mixing process due to the presence of zeros in the wave function becomes effective [26,40]), we expect the relaxation to equilibrium to be accompanied by chaotic rather than stochastic dynamics, as one has in Nelson dynamics (although Nelson dynamics with small but non-zero Brownian motion is hard to distinguish from dBB dynamics, it has the advantage that relaxation is guaranteed to occur, even in the absence of coarse graining and/or mixing).
Ultimately, experiments ought to indicate whether it is relevant, with respect to droplet phenomenology, to formalize the dynamical influence of noise à la Nelson (and/or dBB) as we did in the present paper. We have formulated several proposals in this sense in Sections 6.1 and 8. As emphasized throughout the paper, however, our models should be seen as a first step in the direction of a dynamical model, which remains to be formulated, combining Nelson's stochastic dynamics (and/or dBB dynamics) and memory effects. We think that the results of Section 8 show that this is a promising program for future research.
Finally, it is worth recalling some of the problems that arose when first de Broglie and then Bohm and Nelson developed their theories aimed at deriving quantum dynamics (statistics) as an emergent property, i.e., resulting from an underlying "hidden" dynamics.
The most severe problem is undoubtedly non-locality, which was recognized by Bohm [11,12] to be an irreducible feature of dBB dynamics (see also [19,63] for similar conclusions concerning Nelson-type dynamics). Today, under the influence of the work of John Bell [64] and his followers, it is widely recognized that quantum theory is irreducibly non-local, which makes it particularly difficult to mimic using classical models. Note that entanglement and non-locality (as well as decoherence, which is the corollary of entanglement [65]) only appear if we consider more than one particle at a time, which explains why we did not address these fundamental questions in the core of the paper, where a single droplet is described. It would be interesting to enlarge our model such that the presence of the environment can be taken into account. This would require incorporating the description of open quantum systems, for which a generalization of Bohmian dynamics has been developed in the past [66][67][68][69], but obviously this is beyond the scope of the present paper.
Another problem concerns the fact that the pilot wave is a complex function. This poses still unresolved problems in the case of Nelson dynamics because Nelson's diffusion process does not make it possible [70,71] to fix the phase of the wave function unequivocally (see [72] for an interesting proposal involving a multivalued wave function, also based on Zitterbewegung). In our approach, which is mainly of quantum inspiration, complex wave functions and imaginary phases appear spontaneously, but if we wish to scrutinize the link with the empirically observed modes at the surface of oil baths [9,52,56,61], it will be important to interpret the exact meaning of this complex phase. In the framework of his double solution program [73,74] de Broglie, and others, showed how to derive the Schrödinger equation from a Klein-Gordon equation in the non-relativistic limit. This is only possible provided the real wave bounces at an extremely high frequency (of the order of mc 2 /h). A similar approach has been proposed in the context of droplets phenomenology in [75], where a complex Schrödinger equation is derived from the Klein-Gordon equation along these lines. Although such (interesting and promising) alternative studies of droplets solve the problem of the appearance of a complex phase in a classical context, it is worth noting that the phenomenological results outlined in Section 8, concerning the quantization of droplet orbits in the case of a harmonic potential [9,52], cannot be explained simply in terms of excited modes of the oil bath, because in these experiments only the droplet undergoes the harmonic potential, the oil bath being electromagnetically neutral. This difficulty actually concerns any classical model in which droplet dynamics is formulated in terms of classical modes of the bath only.
To conclude, in our view, the programs that aim at simulating droplet dynamics with quantum tools or at describing the emergence of quantum dynamics based on droplet dynamics, are still largely incomplete and raise challenging fundamental questions. This Pandora's box is now open and it will not be closed any time soon, but this is not something to be feared as it offers new and stimulating perspectives for future research in the field.

Appendix A. Numerical Simulations
Firstly, we discuss the case of the dBB dynamics. It is assumed that we have an analytical solution of the Schrödinger equation Ψ(t, x). We want to compute the evolution of a given initial non-equilibrium density P (t i , x) up to a final time t f and for intermediate time events (we denote all these events by t k , with t 0 = t i and t f = t K ). In particular, we are interested in the coarse-grained non-equilibrium density P (x, t k ) = 1 which is defined in Equation (23). Numerically, we replace that integral by a discrete sum over a finite set of points x l , which are uniformly distributed over the CG cells. In order to obtain the value of each P (x l , t k ) we use the Liouville relation P (x l , t k ) |Ψ(x l , t k )| 2 = where x l i is the position of the particle, which, when evolved according to Equation (11) from t i up to t k , gives x l .
In order to obtain x l i for each x l , we consider the time-reversed evolution with wave-function Ψ * (−t, x) and initial condition x l at time −t k . The position x l , if time evolved from −t k up to −t i according to Equation (11), will give the position x l i . As there is usually no analytical solution of Equation (11), we use a Runge-Kutta (RK) algorithm [76] to obtain a numerical estimate of the position x l i . In order to know if we can trust the result of the RK algorithm, we perform two realizations of the algorithm with different choices of a so-called tolerance parameter (the smaller the value of that tolerance parameter, the more precise the computation), say γ and 10 −1 γ. If the distance between the two positions is less than some chosen value δ, the result of the last iteration of the RK algorithm is trusted. Otherwise, we perform another iteration with 10 −2 γ and we compare it to the previous realization of the RK algorithm. We repeat the procedure until the constraint on the distance between the two successive results of the RK algorithm is satisfied, or until we reach some minimal value of the tolerance parameter. In that case, the position x l is considered as a bad position and it is discarded from the numerical integration of Equation (23). This method was used in [26].
That is one possible method but we could also adopt a more brute-force method: Randomly generate a set of N initial positions according to P (t i , x) and let them evolve according to an Euler algorithm (that is, we divide the time-interval in small time-steps of length ∆t and we increment the position by v(t)∆t at each time step). We record the positions of the N particles for each value of t k , we count the number of particles in each CG cell for each time t k (say n CG ), and we divide n CG by N in order to define P (x, t k ). The first method turns out to be more efficient in the case of the dBB dynamics but it is not applicable in the presence of stochastic terms.
In the case of Nelson dynamics, we used the Euler-Maruyama method for stochastic processes to approximate the solution of the Ito equation, Equation (12). In the same way as Euler's method, the time T is divided into N small discrete time steps ∆t. For each time t i , we generated a random variable normally distributed ∆W i = √ ∆t N (0, 1). The integration scheme has the form We invite the reader interested in the details to consult [77]. The remaining question is how to choose the time step ∆t so that one can trust the result of the numerical simulations. One way to do this is the following. We know that the Born distribution remains invariant under Nelson's dynamics (equivariance). We therefore start with some value for ∆t and decrease it until the result of the numerical simulation confirms this theoretical prediction. We then perform the numerical simulation for the non-equilibrium distribution with the value of ∆t thus obtained.