Timescales for dynamical relaxation to the Born rule

We illustrate through explicit numerical calculations how the Born-rule probability densities of non-relativistic quantum mechanics emerge naturally from the particle dynamics of de Broglie-Bohm pilot-wave theory. The time evolution of a particle distribution initially not equal to the absolute square of the wave function is calculated for a particle in a two-dimensional infinite potential square well. Under the de Broglie-Bohm ontology, the box contains an objectively-existing 'pilot wave' which guides the electron trajectory, and this is represented mathematically by a Schroedinger wave function composed of a finite out-of-phase superposition of M energy eigenstates (with M ranging from 4 to 64). The electron density distributions are found to evolve naturally into the Born-rule ones and stay there; in analogy with the classical case this represents a decay to 'quantum equilibrium'. The proximity to equilibrium is characterized by the coarse-grained subquantum H-function which is found to decrease roughly exponentially towards zero over the course of time. The timescale tau for this relaxation is calculated for various values of M and the coarse-graining length epsilon. Its dependence on M is found to disagree with an earlier theoretical prediction. A power law - tau inversely proportional to M - is found to be fairly robust for all coarse-graining lengths and, although a weak dependence of tau on epsilon is observed, it does not appear to follow any straightforward scaling. A theoretical analysis is presented to explain these results. This improvement in our understanding of timescales for relaxation to quantum equilibrium is likely to be of use in the development of models of relaxation in the early universe, with a view to constraining possible violations of the Born rule in inflationary cosmology.


I. INTRODUCTION
The Born rule is the fundamental connection between the mathematical formalism of quantum theory and the results of experiments. It states that if an observable corresponding to a Hermitian operatorÂ is measured in a system with pure quantum state |Ψ , the probability of an eigenvalue λ i will equal Ψ|P i |Ψ , wherê P i is the projection onto the eigenspace ofÂ corresponding to λ i . For the case of measurement of the position x of -say -an electron in a box, the probability density at time t for finding the electron at x is ρ(x, t) = |Ψ(x, t)| 2 . The Born rule is normally presented as a postulate, though attempts to derive it from more fundamental principles have a long history. There has, for example, been much recent work on deriving the Born rule within the framework of the many-worlds interpretation of quantum mechanics, but such derivations remain controversial [1]. According to a recently-published (2008) encyclopedia of quantum mechanics "the conclusion seems to be that no generally accepted derivation of the Born rule has been given to date, but this does not imply that such a derivation is impossible in principle" [2].
Born's 1926 paper [3], and Heisenberg's introduction of the uncertainty relations the following year [4], were instrumental in popularizing the idea that Nature at the quantum level is fundamentally probabilistic. The idea that Ψ provides a complete description of a single electron strongly suggests that the probabilistic interpretation of |Ψ| 2 expresses an irreducible uncertainty in electron behaviour that is intrinsic in Nature.
It is somewhat ironic therefore -and unknown to most physicists -that Born's rule emerges quite naturally out of the dynamics of a deterministic process that was first outlined by de Broglie in 1927 [5]. The process in question can be described by a theory commonly referred to as the de Broglie-Bohm 'pilot-wave' formulation of quantum mechanics [6][7][8][9][10][11][12][13]. While the theory has attracted little serious interest ever since it was introduced [5,14], there has been a considerable resurgence of activity in this area over the last fifteen years or so [15]. One of the reasons for this is that, although the theory is completely consistent with the full range of predictive-observational data in quantum mechanics, it also permits violations of the Born rule and, at least in principle, this leads to the possibility of new physics and of experimentally testable consequences [16][17][18][19].
The reason that de Broglie-Bohm theory can get away with such an apparently absurd contradiction of one of the basic postulates of quantum theory is that it assumes orthodox quantum mechanics is incomplete, as Einstein always insisted. It supposes that electrons, for example, are real 'particles' with continuous trajectories and that the Schrödinger wave function represents an objectively-existing 'pilot wave' which turns out to influence the motion of the particles. Since the particle density ρ and the square of the pilot wave are logically distinct entities, they can no longer be postulated to be equal to each other. Rather, their identity should be seen as dynamically generated in the same sense that one usually regards thermal equilibrium as arising from a process of relaxation based on some underlying dynamics (though with a dynamics on configuration space rather than phase space). Since pilot-wave theory features a different set of basic axioms and conceptual structures, with event-by-event causality and the prospect of making predictions different from orthodox QM, it is better to think of it as a different theory, rather than a mere 'interpretation' of quantum mechanics.
In the pilot-wave formulation then, quantum mechanics emerges as the statistical mechanics of the underlying deterministic theory. If the particle distribution obeys the Born rule ρ = |Ψ| 2 , the system is said to be in 'quantum equilibrium'. One finds in general that: 1. Non-equilibrium systems naturally tend to become Born-distributed over the course of time, on a coarsegrained level, provided the initial conditions have no fine-grained microstructure [16,[20][21][22]. (The latter restriction is similar to that required in classical statistical mechanics. An assumption about initial conditions is of course required in any time-reversal invariant theory in order to demonstrate relaxation [21,22].) 2. Once in quantum equilibrium a system will remain in equilibrium thereafter, as was originally noted by de Broglie [6] (this property is sometimes referred to as 'equivariance').
These two observations -along with a description of how these statements about the objective makeup of the system might be translated into statements about measurement -can be said to 'explain' or derive the Born rule. Given the common general viewpoint referred to in the first paragraph, many physicists might consider this surprising.
In this work, we present a numerical analysis of the timescale for the relaxation of non-equilibrium distributions of particles to Born-rule quantum equilibrium using pilot-wave dynamics; the approach to equilibrium is monitored by computing the coarse-grained subquantum H -function (see section I B and Ref. 16). The results we obtain are for a particle in a 2D infinite potential square well where the wave function is a finite superposition of M eigenfunctions (where, depending on the choice of initial state, M ranges from 4 to 64).
The initial particle distribution is deliberately chosen to be 'out of equilibrium' by giving it the same form as the absolute square of the ground-state wave function, that is, ρ = 4/π 2 sin 2 x sin 2 y. This system -with fixed M -has been studied before in this context by Valentini and Westman [22] and by Colin and Struyve [23], but here we go further. Our recent development of a new and much faster computer code [24] allows us to study systems with many more modes. The timescale τ for relaxation is studied as a function of the number of modes M (and, in consequence, of the number of nodal points in the wave function) and as a function of the coarse-graining length ε. The dependence of the relaxation timescale on these two quantities is compared to theoretical predictions.
It is intended that calculations such as these will provide a next step towards a detailed understanding of relaxation to quantum equilibrium in the early universe, with a view to constraining possible non-equilibrium effects in cosmology. In Ref. 19 it was shown, in the context of inflationary cosmology, that corrections to the Born rule in the early universe would in general have potentially observable consequences for the cosmic microwave background (CMB). This is because, according to inflationary theory, the primordial perturbations that are currently imprinted on the CMB were generated at early times by quantum vacuum fluctuations whose spectrum is conventionally determined by the Born rule. To make detailed predictions for possible anomalies in the CMB, however, requires a precise understanding of how fast relaxation would occur in, for example, a pre-inflationary era (as discussed in section IV-A of Ref. 19). It may be hoped that numerical studies, such as those reported in this paper, will reveal how the relaxation timescale depends on general features of the quantum state such as the number M of modes in a superposition. The results could then be applied in future work to specific cosmological models.

A. Pilot-wave dynamics
The basic ideas of de Broglie-Bohm pilot-wave theory may be simply understood in a non-relativistic context [25]. It is a non-local hidden-variables theory, that is, the theory contains some variables that distinguish the individual members of an ensemble that in orthodox QM would be considered identical since they all have the same wave function. These variables are supposed to be ultimately responsible for the apparently random nature of -for example -position measurements on the system. If, as required by some interpretations, one were to suppose both that a complete description of the system is afforded by Ψ and that Ψ has an objective, physical existence, one might conclude from the results of measurements that Nature is intrinsically probabilistic or random. In pilot-wave theory, by contrast, one supplements the wave function description with 'hidden variables' by postulating the existence of particles with definite positions, in addition to the wave.
These particles then follow deterministic trajectories (the nature of which can be deduced) and the observed randomness is then understood to be a consequence merely of our ignorance of the initial conditions, that is, the starting positions of the particles.
How does an individual quantum system evolve in time? The pilot wave evolves at all times according to the usual time-dependent Schrödinger equation As normally understood the evolving quantum system behaves like a 'probability fluid' of density |Ψ| 2 = ΨΨ * with an associated time-dependent quantum probability current, defined in the usual manner as j = m Im(Ψ * ∇Ψ). In pilot-wave theory, the particles have a continuous objective existence, with trajectories that follow the streamlines of the current. Thus their velocity is given by the current divided by the density, that Using the complex polar form of the wave function Ψ = |Ψ| exp[iS/ ], we can recover the (locally defined) phase S(x 1 , . . . , x N , t) of the wave by the expression S = Im ln Ψ. The de Broglie guidance equation for the trajectories x i (t) may then be written as If, for an ensemble of particles with the same wave function, the initial positions have a Born-rule distribution, then (by construction) the law of motion of Eqn. 1 implies that the particle positions will have a Born-rule distribution at all times.
If desired, one may take the first time derivative to write the equation of motion in second-order form, where the quantum potential In this approach, the system acts as if there were a 'quantum force' −∇ i Q acting on the particles in addition to the classical force −∇ i V . This second-order approach with a law of motion given by Eqn. 2 was proposed by Bohm in 1952. It may be referred to as 'Bohm's dynamics' in order to distinguish it from 'de Broglie's dynamics' based on Eqn. 1 (which was proposed by de Broglie in 1927). For de Broglie, p = ∇S is the law of motion; for Bohm -at least the Bohm who wrote the 1952 papers -it is an initial condition which can be dispensed with (clearly if we integrate the second-order formula we only recover de Broglie's equation up to some constant and this must be fixed for each trajectory by some boundary condition, such as that implied by de Broglie's equation for some time t 0 ). Thus, in principle, Bohm's dynamics encompasses what one might call 'extended nonequilibrium' where p = ∇S in addition to ρ = |Ψ| 2 . Recent work [26] suggests that this 'extended nonequilibrium' is unstable and does not relax in general; if this is correct then it may be argued that Bohm's second-order dynamics is untenable as a fundamental theory as there would be no reason to expect equilibrium in the universe today, and that de Broglie's dynamics is in fact the fundamental formulation of pilot-wave theory.
Some additional relevant observations: (1) The form of the guidance equation may be altered, while retaining consistency with the Born-rule distribution. This can be achieved by adding a divergence-free term (divided by |Ψ| 2 ) to the right-hand side. Such alternative velocity fields will not be discussed further here but have been studied by, for example, Colin and Struyve [23] and Timko and Vrscay [27]. Note that such alternatives yield an equivalent physics only in the equilibrium state; away from equilibrium, 'subquantum' measurements would allow one to track the trajectories and so distinguish the true velocity field [17].
(2) Given the wave function for a system, the particle trajectories from any starting point may be calculated using only the initial position of the particle, rather than the position and the momentum. This is because the guidance equation alone gives the particle velocity and consequently the momentum for any initial position.
(3) Particle trajectories tend to be quite erratic, even with simple wave functions that are superpositions of just a few energy eigenfunctions. Fig. 1 illustrates the divergence of neighbouring particle trajectories by showing the paths of two particles with almost identical initial positions, propagating according to pilot-wave dynamics.
How do numerical simulations demonstrating the Born rule for the actual particle positions translate into statements about 'measurement' ? Ideal measurements of position in pilot-wave theory are usually correct measurements (they reveal the pre-existing position of the particle -see Ref. [9], p.351) and so the Born rule in position space follows immediately if the particles really are distributed that way. For other kinds of measurements, a clear derivation of the Born rule may be found in Section 8.3.5 of Holland's textbook [9] (noting that Holland assumes the |Ψ| 2 distribution of actual particle positions as a postulate -see Ref. 9, p.67). The key point is that, in a theory of particles, experimental observations may be reduced to particle positions (dots on screens, apparatus pointer positions, etc.) -where laboratory apparatus is treated as just another system made of particles. As long as the Born rule holds for the joint distribution of positions of all the particles involved (including the particles making up the equipment), then the marginal probability distribution for, say, pointer positions (obtained by integrating out the other degrees of freedom) will necessarily match the predictions of quantum mechanics. In such a case, the distribution of macroscopically-recorded outcomes will be the same as in quantum theory. The particles have been propagated for time 4π in both cases -note the rapid divergence. 1(a) also shows a good example of particle motion near to a node. The trajectory is seen to spiral around a moving nodal point before departing from the vicinity of the node (similar behaviour is reported in Ref. 22). This behaviour seems to be a major driving force behind relaxation.

B. Quantum equilibrium
To demonstrate equivalence to quantum mechanics, it is usually simply assumed that the actual distribution of particle positions is already supplied to us obeying the Born rule ρ = |Ψ| 2 . In the approach taken here, where we try to demonstrate why this is so, the Born-rule distribution is considered to be a special case and the particles are said to be in quantum equilibrium when in this state. The dynamics described in section I A can just as well be used to describe the evolution of non-equilibrium systems, whereas standard formulations cannot. In general in such studies, the probability density is found to approach the Born rule distribution over time; it is said to relax to equilibrium [16,22]. This relaxation is a consequence of the deterministic motion of the particles and is not an intrinsically stochastic process (further insight into relaxation has been obtained by Bennett using techniques from Lagrangian fluid dynamics [28]). Fig. 2 shows the results of a numerical simulation of this relaxation process. It can be clearly seen that the particle distribution ρ rapidly comes to resemble the (periodically repeating) time-dependent |Ψ| 2 . The example chosen -a superposition of sixteen modes, for a particle moving in two spatial dimensions -is identical to that studied by Valentini and Westman [22]. The results obtained match theirs, thereby providing an important confirmation of the previous results, with an independently written and implemented numerical code [29]. Computational model of the relaxation of an initially non-equilibrium distribution, ρ = 4 π 2 sin 2 x sin 2 y, evolving according to pilot-wave dynamics. The wave function is a superposition of the first sixteen eigenstates for a particle in a 2D infinite potential square well. The simulation was run for one period of the wavefunction, or 4π in these units. Even after such a short period, significant relaxation towards equilibrium can be observed. (These results provide an independent confirmation of those first obtained by Valentini and Westman [22].) If we are to have any chance of observing new physics associated with quantum nonequilibrium states [16][17][18][19][20][21]30], it is important to have a handle on the timescale of this relaxation. To quantify the proximity of a distribution to equilibrium, we may use an analogue of the classical H -function [16,22]. This 'subquantum H -function' is defined as where the (weighted) volume element dΣ = |Ψ| 2 dx.
This quantity will be zero if and only if ρ = |Ψ| 2 everywhere, and will be positive otherwise [31], making it a useful measure of proximity to equilibrium. Clearly, H is simply the negative of the relative entropy of ρ with respect to |Ψ| 2 .
A feature of the quantities in this definition is that both the volume element, dΣ, and the ratio f = ρ |Ψ| 2 are preserved along trajectories. To show this for f , consider the two continuity equations: which follows from the assumption that the actual trajectories follow the velocity field given by Eq. 1, and which follows from the Schrödinger equation.
These two equations can be used to show that the ratio f = ρ |Ψ| 2 obeys: Thus f will be preserved along trajectories. Thus, if the system is initially in quantum equilibrium, with f = 1 everywhere, it will never depart from that state. This can, of course, be seen directly from the fact that ρ and |Ψ| 2 obey identical continuity equations: if ρ and |Ψ| 2 are initially equal, they will necessarily remain equal at all times, since their time evolutions are determined by the same partial differential equation.
For general (non-equilibrium) initial conditions, the exact value of the H -function remains unchanged as the system evolves. However, if a coarse-graining is applied to ρ and |Ψ| 2 , that is, we replace ρ → ρ, |Ψ| 2 → |Ψ| 2 where overbars indicate averaging over small coarse-graining cells, then the coarse-grained H -function can be shown to be non-increasing, on the assumption that the initial state contains no fine-grained microstructure (as in the analogous classical coarse-graining H -theorem) [16]. Furthermore, H will in fact decrease, if the initial velocity field varies with position across the coarse-graining cells [20,21]. The decrease of H represents a relaxation of the system towards equilibrium, and formalizes an analogue of the intuitive idea of Gibbs: an initial non-equilibrium distribution will tend to develop fine-grained microstructure and become closer to equilibrium on a coarse-grained level. Heuristically speaking, this may be thought of in terms of two 'fluids', with densities ρ and |Ψ| 2 , that are 'stirred' by the same velocity field, and thereby tend to become indistinguishable when coarse-graining is applied.
The effects of coarse-graining on the particle density at some randomly-selected time may be seen in Fig. 3. 3(a) shows a snapshot of the fine-grained density ρ on a 1024 × 1024 lattice; 3(b) shows the coarse-grained density derived from averaging the fine-grained density over square cells containing 32 × 32 lattice points (in which case we say the system has a 'coarse-graining length ε = 32'). 3(c) shows the result of a 'smoothed' coarse-graining -using overlapping cellswhich is more suitable for plotting graphs [22]. 3(d) is the same coarse-grained density as in Fig. 3 For the 16-mode case, it was found in Ref. 22 that H decays approximately exponentially as One of us (AV) has presented a theoretical estimate of the relaxation timescale τ obtained by considering the behaviour of the second-time derivative of H at t = 0 (where H posesses a local maximum) [20,21].
As discussed in more detail below, it was shown that in the limit ε → 0, τ scales inversely with ε. Further estimates, or simply dimensional analysis, then suggested the rough formula [21] τ Here ε is the length of the coarse-graining cells and ∆E is the energy spread of the wave function. For reasonable values of ε, this estimate was in rough agreement with the numerical value [22].
However, because H has a local maximum at t = 0, the estimate on the right-hand side of Eqn. 9 -obtained from the second time-derivative of H at t = 0 -can only define a timescale that is valid close to t = 0. As we shall see, it cannot properly represent the timescale τ associated with the subsequent (approximately) exponential decay. The results of this paper in fact show a scaling of τ with ∆E that disagrees with Eqn. 9 -but which is in agreement with an improved estimate discussed below.
The question of how τ varies with the number of modes was not investigated by Valentini and Westman [22] owing to computational difficulties. That gap is filled in this paper.

II. NUMERICAL SIMULATIONS
In this work we compute the dependence of the relaxation time τ on the coarse-graining length ε and energy spread ∆E through explicit numerical simulations. An initially non-equilibrium probability density in a 2D infinite potential square well is evolved according to pilot-wave dynamics, using a wave function consisting of an out-of-phase superposition of the first M energy eigenstates (normal modes). For this choice of wave function, taking all the modes to have equal weight, we have ∆E ∼ M 2 so we look for a dependence of the On the basis of Eqn. 9, for example, we would expect p = −3 and q = −1.
To study this system (and potentially others, since it is designed to be easily extendible) we have written a new computer code named 'LOUIS' [24] which uses pilot-wave dynamics to calculate particle trajectories.
Given an initial probability density and wave function, LOUIS is able to use these trajectories to compute the probability density function and coarse-grained subquantum H -function at later times. It is a ground-up reimplementation of the code used in Refs. 22 and 23 and is up to two orders of magnitude faster with significantly more capabilities. Currently it can treat infinite potential square wells in one, two, or three dimensions using a finite superposition of eigenstates to represent the wave function. The relative weights and phases of the eigenstates may be specified in the input or chosen randomly (but reproducibly, using preset seeds). The scale of the coarse-graining may also be set manually, outputting results for multiple coarse-graining lengths in a single run of the program.
Since the results of interest involve calculation of the subquantum H -function, which in turn involves a numerical integration over the area of the potential well, the quantities in the integrand, ρ and |Ψ| 2 , must be evaluated on a regular lattice. In all calculations presented here, a 1024 × 1024 lattice is used covering a square two-dimensional cell of length π. The 'coarse-graining length' ε refers to the number of lattice points along one side of a coarse-graining cell.

A. Details of the algorithm
The LOUIS code uses de Broglie-Bohm trajectories to calculate how the particle probability density evolves from a given initial density. At each of a sequence of requested times, it evaluates the particle density and wave function at all points on the fine-grained lattice, and then applies coarse-graining on the requested scales.
The coarse-grained H -function is calculated from these data at each timestep, and output files containing How do we calculate the density at a later time? We have seen that the ratio ρ |Ψ| 2 is preserved along trajectories, implying that the density at position x and time t may be calculated from where the positions x 0 and x are points on the same trajectory, at times 0 and t respectively. The value of |Ψ| 2 can be calculated analytically at all positions and times, and ρ(x 0 , 0) is a known function, therefore ρ(x, t) may be calculated directly once we know the trajectory endpoint x. This is the crucial relation used to calculate probability density functions from trajectories.
In fact, certain practicalities require real calculations to be performed in a slightly different manner. The subquantum H -function is evaluated through numerical integration over the 2D box from a set of values of ρ and |Ψ| 2 calculated at discrete points. Since accurate and efficient quadrature algorithms in few dimensions generally require the points to be sampled uniformly across the region, LOUIS starts with a uniform lattice at time t and exploits the time-reversibility of the dynamics to calculate particle trajectories backwards in time to t = 0. This ensures uniform sampling of ρ(x, t) at t when the quadrature is to be performed. This has the unfortunate consequence that if ρ is required at a later time t > t this 'backtracking' has to be done all the way to t = 0 again: the data calculated at time t cannot be used again.
The rate-limiting step of the LOUIS program is the numerical integration of the de Broglie guidance equation v(x, t) = ∇S(x, t) (in atomic units) to compute the particle trajectories x(t). One may use a variety of standard algorithms; an excellent choice for these purposes is Runge-Kutta-Fehlberg [32]. Currently, the Schrödinger equation is not integrated numerically to compute the time-development of the wave function; instead, only finite superpositions of stationary states are used so the wave function can be evaluated exactly for any t.
The velocity of the particle at any point may be computed from Im ∇Ψ(x,t) Ψ(x,t) where the M -mode wave function is given by Here E mn are the energy eigenvalues 1 2 (m 2 + n 2 ), the θ mn are the (randomly chosen) initial phases, m, n = 1, 2, · · · , √ M are positive integers, and (for convenience) M has an integer square root.
As with all such algorithms, Runge-Kutta-Fehlberg basically involves adding small increments to a functionhere x(t) -where the increments are given by derivatives ( dx dt = v = ∇S = Im ∇Ψ Ψ ) multiplied by variable step sizes (here, a timestep ∆t). In order to increase the accuracy, a tolerance is set for the maximum error on each step (the step tolerance); if the error is greater than this, a smaller timestep is used (subject to appropriate underflow checks). When the integration has been performed along the entire trajectory between the required initial and final times, the whole trajectory is recomputed with the step tolerance decreased by a factor of 10. If the two final positions agree within a certain tolerance (the trajectory tolerance), then the trajectory is kept. If not, the process is repeated with smaller and smaller step sizes until the trajectories converge, or until the step tolerance reaches a certain minimum value where the calculation will take too much time and the trajectory is flagged as failed. Failed trajectories are not used in the subsequent computation of the density. The proportion of failed trajectories rarely exceeds 1 in 1000 and their contribution to the overall error is negligible. In general, they are trajectories that come too close to a wave function node (where the velocity field diverges).
Computational cost is the main limiting factor. The calculation of relaxation timescales is very computationally intensive, requiring many long, high-precision numerical integrations. Since the particle wave function and its gradient must be evaluated at each step in the integration of each trajectory, the complexity of the wave function contributes significantly to the time taken to perform the calculation; runs with larger numbers of modes in the superposition are considerably more expensive. For example, a typical calculation on an elderly cluster of sixteen processors (seventeen evaluations of H between time 0 and time 4π) took 9.6 CPU-hours with a 9-mode wave function. A comparable calculation with 36 modes took 542 CPU-hours.

III. RESULTS AND DISCUSSION
We begin by verifying that exponential decay is an appropriate model for the evolution of H under these conditions as was assumed in the definition of the relaxation timescale. Fig. 4 shows plots of ln H vs. t for different coarse-graining lengths ε, different numbers of modes M in the superposition, and different sets of initial phases θ mn in the wave function (in order to obtain reproducible results the phases are fixed by a single 'preset' parameter in the LOUIS input file, which controls the seed for the random generation of a set of phases). In all cases we find a good straight-line fit to the data, validating the assumption of exponential decay in H over this range of conditions. This was not unexpected as exponential decay was previously demonstrated by Valentini and Westman [22], though only for the case of M = 16 modes with a fixed coarse-graining length.
A. Relaxation time τ as a function of the number of modes M Fig. 5 shows plots of ln τ against ln M for various coarse-graining lengths ε where now logarithmic axes are used in order to search for a power law relationship of the form τ ∝ M p . Error bars were calculated by running LOUIS six times for each (ε, M ) pair with different initial phases θ mn in each run; the mean of these was taken to be the best estimate for the timescale with the standard deviation taken to be the error bar.
In the case of 4-mode simulations, the spread of values of τ was so large that the error bar ∆τ could not be displayed on a logarithmic plot (hence the arrow at the base of the corresponding error bars in Fig. 5).
Some representative values of τ were (980 ± 1600) for ε = 4 and (100 ± 110) for ε = 32. This large spread in the timescale for low M can be understood by considering the role of wave function nodes in the relaxation.
The rapidly-varying velocity field in the vicinity of nodes is believed to be a significant driving force for this process [22], and so the initial positions of the nodes are likely to affect the timescale (these initial positions are moved around when modifying the set of initial phases θ mn ). The important point is that in larger The horizontal axis on the first graph in 4(b) is scaled differently to the others, to better fit the range of data.
superpositions with larger M there are many nodes and the exact change in positions of nodes will have less effect because the average distribution will be similar. With a small superposition, maybe only containing one node, the initial position of this node will have a much larger effect on the subsequent relaxation, and so the six runs with different initial phases will tend to produce very different results.
A full animation of the relaxation process can be seen online at Ref. [33]. The swirling vortices in the density surrounding the moving wave function nodes are quite striking, and one can obtain more of a visual sense of why the presence of nodes increases the chaotic nature of the trajectories.
In Fig. 6, therefore, we choose to plot the same data as in Fig. 5  The approximation used in this analysis is that, ∆E ∼ M 2 , which only applies for larger M .
The best-fitting power law is shown on the graph in each case, and they are also summarised in Table I. Clearly these results do not support the theoretical prediction (described as a 'crude estimate' in Ref. [22]) that the relaxation time should be proportional to M −3 . Instead they very strongly suggest a relationship of the form τ ∝ M −1 , to within the estimated error. This suggests that some of the approximations made in obtaining the prediction in Eqn. 9 are invalid.
The relevant arguments used to obtain the apparently incorrect estimate of the scaling are set out in Ref. [21].
They begin by defining a relaxation timescale τ in terms of the rate of decrease of H near t = 0 via the following formula: 1 The second derivative is used rather than the first derivative since dH dt 0 = 0 (as may be shown analytically [20]. The H-curve necessarily has a local maximum at t = 0. This property might, at first sight, seem to be incompatible with the observed exponential decay. But in fact, it must be the case that the exponential decay sets in soon after t = 0, and that in the limit t → 0 the decay is not exponential. The timescale of Eqn. 13 then applies only to this very short period immediately after t = 0. In other words, while Eqn. 13 may well be a good estimate for the relaxation timescale in the limit t → 0, it cannot accurately estimate the time constant in the exponential tail -where the bulk of the relaxation takes place -and is therefore of little practical relevance. Another potential conflict in the derivation of Eq. 9 is the requirement for a logarithmic scale, as the lower bound is less than zero. There is reason to believe (see section III) that these points may be excluded, and with such a large error their weight in fitting would be very small, so Fig. 6 shows the same results without these points.  coarse-graining length ε so small that the velocity field varies little over the length of a coarse-graining cell.
The derivation works by considering the dependence of τ on ε in the limit where ε → 0, then applying dimensional analysis to find the other dependencies. The inverse dependence of τ (as defined by Eqn. 13) on ε must hold in this limit, since we have an analytic proof of it. However, unfortunately, the limit seems too restrictive to be of practical use. In the cases studied here the wave function apparently varies rapidly enough on that scale -particularly with larger numbers of modes in the superposition -to ensure that we are not working in this limit at all. For example, Fig. 7 demonstrates the effect of coarse-graining (ε = 32) on one of our 64-mode wave functions. This shows the magnitude of Ψ, rather than the velocity field, but the length scale over which the latter varies should be even smaller than the length scale over which the former varies. We may thus conclude that the velocity field is likely to vary significantly over the length of one coarse-graining cell, at least under some of the conditions studied here. Taking all these facts into consideration, it is not surprising that the results do not confirm Eqn. 9, whose domain of validity is probably simply too narrow to be of practical use.
We now provide a theoretical justification for the observed relation τ ∝ M −1 . Let q r (r = 1, 2) denote x or y. We shall proceed by considering an upper bound on the (equilibrium) mean displacement δq r of particles over an arbitrary time interval [t i , t f ]. A relaxation time may then be defined, in the case of the infinite potential square well, by the condition that relaxation will occur over timescales τ such that the said upper bound becomes of the order of the width L of the potential well. As we shall now show, the timescale will then be τ ∼ L m 2Ē (14) whereĒ is the mean energy, andĒ ∝ M 2 .
There now follows a derivation of the upper bound on the mean displacement δq r . The derivation is based on Ref. 34, but with some differences that are highlighted below.
First, note that the final displacement δq r (t f ) has modulus whereq r is the component of the de Broglie-Bohm velocity in the q r direction.
Using the fact that, for any x, x ≤ x 2 , we have Now note that, starting from the guidance equation: whereπ r is the momentum operator conjugate to q r , and π r denotes the usual quantum expectation value for the operatorπ r . The last equality follows from the relations and Since (∂ |Ψ| /∂q r ) 2 ≥ 0, we then have and so We also have π2 where W r denotes the x-or y-part of the mean Hamiltonian, with W r ∝ M 2 . Hence, Since W r is time-independent, and setting t i = 0 and t f = t we have Setting the right-hand side to be of order L, and noting thatĒ ≈ W r , then indeed yields the relaxation time in Eqn. 14, whose inverse scaling with M is in agreement with the numerical results presented above.
As mentioned, this derivation is based on that in Ref. 34 but differs in some respects. The purpose of the analysis in Ref. 34 was to derive a condition for the suppression of relaxation in expanding space (here we are only concerned with static space) and the condition for relaxation was that the mean displacement δq rfor field degrees of freedom in Fourier space -should be comparable to (or greater than) the quantum spread ∆q r . In the analysis above, the only degree of freedom considered is the spatial displacement of a particle in the potential well, the constraints of which slightly change the condition for relaxation. Regardless of the spread in the wave function the particle cannot move beyond the confines of the well, so the condition used for relaxation is that the mean displacement of a particle is comparable to (or greater than) the size of the potential well.   Fig. 8.
The errors are estimated in the same way as the results in Table I. These values are not compatible with Eqn. 9, τ ∝ ε −1 , nor do they appear to be compatible with a consistent power law.

B. Relaxation time as a function of coarse-graining length
In Fig. 8, the relaxation time τ is plotted against coarse-graining length ε for various M , again with logarithmic axes and with the data for M = 4 excluded for the same reason as before. The best-fit power law is shown on the graphs, and the results are summarized in Table II. It is evident from these data that, although the relaxation time evidently decreases with increasing ε, the data do not support Eqn. 9 nor are they particularly suggestive of a constant power law. A weak dependence of order τ ∼ ε − 1 / 4 is observed.
As was the case with the M -dependence of τ , the dependence on ε is not in fact expected to take the form in Eqn. 9, for two reasons. First, the derivation of Eqn. 9 is based on a definition of 'relaxation time' that applies only when t → 0; it does not apply to the exponential tail of H, where most of the relaxation takes place. Second, the velocity field can vary significantly over a coarse-graining cell, contrary to the assumption made in the derivation of Eqn. 9. Indeed the graphs in Fig. 8 appear to show a systematic concave character rather than a straight line, apparently suggesting a small systematic deviation from a power law model. The concave curvature is at least consistent with a power law τ ∝ ε −1 in the limit ε → 0 (which must hold as we have noted), since at some point to the left of the data in any of the figures the gradient ought to approach (or pass through) −1. A systematic study at smaller coarse-graining lengths would however be rather difficult in terms of CPU time because of the need for a significantly finer lattice and much greater overall number of lattice points.
The lack of a dependence on ε in Eqn. 14 is not surprising since the analysis in section III A and [34] uses a different definition of the timescale. This definition does not consider theH-function nor is there any necessity to even mention coarse-graining. The weak dependence observed in the numerical simulations should be interpreted as an effect outside the scope of this prediction rather than one in conflict with it.

IV. CONCLUSIONS
The numerical simulations performed in this work demonstrate clearly and unequivocally the tendency for Born-rule distributions to arise spontaneously as a consequence of ordinary pilot-wave dynamics, even for a system as simple as the electron in a two-dimensional potential well. Contrary to popular belief therefore, prediction, Eqn. 9, but agrees with an alternative theoretical analysis presented here in section III A. As we have discussed, the assumptions made in deriving Eqn. 9 were probably too restrictive for it to be of practical use. In particular, the defined timescale is relevant only close to t = 0, and does not apply to the exponential tail of the H function, where most of the relaxation takes place.
Our simulations reveal no well-defined scaling for the relaxation time as a function of coarse-graining length τ (ε), other than a possible weak-dependence of the order of τ ∼ ε − 1 / 4 . This also differs from Eqn. 9 and this is probably simply because the coarse-graining cells are too large for the derivation of Eqn. 9 to be valid. It is possible that by decreasing ε a behaviour conforming better to Eqn. 9( τ ∝ ε −1 ) could be observed, and it would be interesting to see at what length scale this begins to emerge.
Physically speaking our results suggest very short relaxation times with a range of values observed for τ between about 1 and 1000. Using natural units c = = 1 and an electron with mass m = m e = 1 this corresponds to relaxation times of the order of 10 −21 -10 −18 s. This is consistent with our current understanding of quantum mechanics and modern experimental investigations in which no deviation from quantum equilibrium is observed. If the initial state of the universe corresponded to a non-equilibrium state, one might then assume that almost all deviations from the Born rule will have been quickly washed out.
However, it may be (see Refs. 19 and 34 and works cited therein) that relic non-equilibrium from the early universe could be observed today, either directly or by its imprint on the cosmic microwave background. In the future we intend to modify the LOUIS program to deal with more realistic wave functions, multi-particle systems, and expanding spaces with the intention of improving some of the predictions outlined in Refs. 19 and 34. More precise predictions for the effect on the CMB of non-equilibrium in the early universe could lead to experimental tests of the de Broglie-Bohm formulation of quantum mechanics.
How might these results generalize to more complex systems? In the original version of the subquantum H -theorem published by one of us (AV) in 1991 [16], relaxation was considered for a theoretical ensemble of complex N -body systems. Once equilibrium is reached for such an ensemble, it can be (and was) shown that extracting a single particle from each system resulted in single-particle sub-ensembles that obey the Born rule. The original view expressed in Ref. 16 was that, realistically, relaxation would take place efficiently only for many-body systems, and that the Born rule for single particles would be derived by considering how these are extracted from more complex systems. But in fact, in practice, there is an efficient relaxation even in simple two-dimensional one-electron systems, and there appears to be no need to appeal to a complex N -body 'parent system'.
Note that in a strict account of our world, say in the early universe, it could well be that all degrees of freedom are entangled, so that there is no actual ensemble of independent subsystems with the same wave function. However, one can still talk about a theoretical ensemble of universes, each with the same universal wave function, and consider the evolution of its distribution. (One could also consider a mixed ensemble of universes, and apply our discussion to each pure sub-ensemble.) It is sometimes suggested that it is problematic to consider probabilities for the 'whole universe'. And yet, cosmologists are currently testing primordial probabilities experimentally by measuring temperature anisotropies in the cosmic microwave background. By making statistical assumptions about a theoretical 'ensemble of universes', cosmologists are able to test probabilities in the early universe, such as those predicted by quantum field theory for vacuum fluctuations during inflation. (For a detailed discussion of this in a de Broglie-Bohm context, see Ref. [19].) One can question what the ensemble of universes refers to. Is it a subjective probability distribution? Or, is the universe we see in fact a member of a huge and perhaps infinite ensemble, as is the case in theories of eternal inflation? Those are interesting questions, but only tangentially related to the ongoing experimental tests. It is also important to bear in mind that there is much that is not known about cosmology, so the treatment should be kept independent of cosmological details as far as possible. What we do know is that all the particles we see today are, or were, in complex superpositions (whether entangled with other particles or not), and it appears clear from the simulations that such superposition yields rapid relaxation -if it is rapid in two dimensions, one would expect it to be even more rapid for 3N dimensions.