## Abstract

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 Schrödinger 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 time scale *τ* for this relaxation is calculated for various values of *M* and the coarse-graining length *ε*. Its dependence on *M* is found to disagree with an earlier theoretical prediction. A power law, *τ*∝*M*^{−1}, is found to be fairly robust for all coarse-graining lengths and, although a weak dependence of *τ* on *ε* 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 time scales 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.

## 1. 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 |*Ψ*〉, then the probability of an eigenvalue *λ*_{i} will equal , where 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-world interpretation of quantum mechanics, but such derivations remain controversial (Saunders *et al.* 2010). According to a recently published encyclopaedia 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’ (Landsman 2008).

Born's (1926) original paper and Heisenberg's (1927) introduction of the uncertainty relations the following year 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 (Bacciagaluppi & Valentini 2009). The process in question can be described by a theory commonly referred to as the de Broglie–Bohm ‘pilot-wave’ formulation of quantum mechanics (de Broglie 1928; Bohm 1952*a*,*b*; Bohm & Hiley 1993; Holland 1993; Dürr & Teufel 2009; Riggs 2009; Towler 2009; Valentini 2009). While the theory has attracted little serious interest ever since it was introduced—possible reasons why being discussed in Cushing (1994), Bacciagaluppi & Valentini (2009) and Towler (2009)—there has been a considerable resurgence of activity in this area over the last 15 years or so (Towler & Valentini 2010). 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 (Valentini 1991*a*,*b*, 2002*b*, 2007, 2010).

The reason that the 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 on 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 quantum mechanics, 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}, then the system is said to be in ‘quantum equilibrium’. One finds in general that:

— Non-equilibrium systems naturally tend to become Born-distributed over the course of time, on a coarse-grained level, provided the initial conditions have no fine-grained microstructure (Valentini 1991

*a*,*b*, 1992, 2001; Valentini & Westman 2005). 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 (Valentini 2001; Valentini & Westman 2005).— Once in quantum equilibrium a system will remain in equilibrium thereafter, as was originally noted by de Broglie (1928). This property is sometimes referred to as ‘equivariance’.

These two observations—along with a description of how these statements about the objective make-up 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 *time scale* 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 §1*b* and Valentini 1991*a*,*b*). The results we obtain are for a particle in a two-dimensional 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, . This system—with fixed *M*—has been studied before in this context by Valentini & Westman (2005) and by Colin & Struyve (2010), but here we go further. Our recent development of a new and much faster computer code (Towler 2010) allows us to study systems with many more modes. The time scale *τ* 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 time scale on these two quantities is compared with 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 Valentini (2010) 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 Valentini (2010, §4*a*). It may be hoped that numerical studies, such as those reported in this paper, will reveal how the relaxation time scale 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.^{1} It is a non-local hidden-variables theory, that is, the theory contains some variables which distinguish the individual members of an ensemble that in orthodox quantum mechanics 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 . 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 is, by
Using the complex polar form of the wave function , we can recover the (locally defined) phase *S*(**x**_{1},…,**x**_{N},*t*) of the wave by the expression . The *de Broglie guidance equation* for the trajectories **x**_{i}(*t*) may then be written as
1.1If, 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 equation (1.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 the second-order form,
1.2where 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 equation (1.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 equation (1.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 non-equilibrium’ where **p**≠∇*S* in addition to *ρ*≠|*Ψ*|^{2}. Recent work (Colin *et al.* in preparation) suggests that this ‘extended non-equilibrium’ 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:

— 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 & Struyve (2010) and Timko & Vrscay (2009). 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 (Valentini 2002*b*).— 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.

— Particle trajectories tend to be quite erratic, even with simple wave functions that are superpositions of just a few energy eigenfunctions. Figure 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 Holland 1993, 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 Holland (1993, §8.3.5) (noting that Holland assumes the |*Ψ*|^{2} distribution of actual particle positions as a *postulate* (see Holland 1993, 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.

### (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 §1*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 (Valentini 1991*a*; Valentini & Westman 2005). 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 using techniques from Lagrangian fluid dynamics (Bennett 2010)).

Figure 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 16 modes for a particle moving in two spatial dimensions—is identical to that studied by Valentini & Westman (2005). The results obtained match theirs, thereby providing an important confirmation of the previous results, with an independently written and implemented numerical code.^{2}

If we are to have any chance of observing new physics associated with quantum non-equilibrium states (Valentini 1991*a*, 1992, 2001, 2002*a*,*b*, 2007, 2010), then it is important to understand the time scale of this relaxation. To quantify the proximity of a distribution to equilibrium, we may use an analogue of the classical *H*-function (Valentini 1991*a*,*b*; Valentini & Westman 2005). This ‘subquantum *H*-function’ is defined as
1.3

This quantity will be zero if and only if *ρ*=|*Ψ*|^{2} everywhere, and will be positive otherwise,^{3} 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 this definition is that the ratio *f*=(*ρ*/|*Ψ*|^{2}) is preserved along trajectories. To show this for *f*, consider the two continuity equations
1.4which follows from the assumption that the actual trajectories follow the velocity field given by equation (1.1), and
1.5which follows from the Schrödinger equation.

These two equations can be used to show that the ratio *f*=(*ρ*/|*Ψ*|^{2}) obeys
1.6

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 , where the overbar indicates averaging over small coarse-graining cells, then the coarse-grained *H*-function
1.7can be shown to be non-increasing (Valentini 1991*a*), on the assumption that the initial state contains no fine-grained microstructure (as in the analogous classical coarse-graining *H*-theorem). Furthermore, will in fact decrease, if the initial velocity field varies with position across the coarse-graining cells (Valentini 1992, 2001). The decrease of 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 figure 3.

For the 16-mode case, it was found in Valentini & Westman (2005) that decays approximately exponentially as
1.8One of us (A.V.) has presented a theoretical estimate of the relaxation time scale *τ* obtained by considering the behaviour of the second time derivative of at *t*=0, where possesses a local maximum (Valentini 1992, 2001). As discussed in more detail below, it was shown that, in the limit , *τ* scales inversely with *ε*. Further estimates, or simply dimensional analysis, then suggested the rough formula (Valentini 2001)
1.9Here, *ε* 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 (Valentini & Westman 2005).

However, because has a local maximum at *t*=0, the estimate on the right-hand side of equation (1.9)—obtained from the second time derivative of at *t*=0—can only define a time scale that is valid close to *t*=0. As we shall see, it cannot properly represent the time scale *τ* associated with the subsequent (approximately) exponential decay. The results of this paper in fact show a scaling of *τ* with Δ*E* that disagrees with equations (1.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 & Westman (2005) owing to computational difficulties. That gap is filled in this paper.

## 2. 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 two-dimensional 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 form
2.1On the basis of equation (1.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’ (Towler 2010) 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 codes used in Valentini & Westman (2005) and Colin & Struyve (2010) 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 time step, and output files containing (*t*,*H*) pairs and pairs are written to disk. The program calculates a straight-line fit using linear regression of the *t* versus data; assuming exponential scaling, the gradient of this is the decay constant or relaxation time *τ*.

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
2.2where 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 two-dimensional 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 the Runge–Kutta–Fehlberg algorithm (Press *et al.* 1992). 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
2.3Here, *E*_{mn} are the energy eigenvalues , the *θ*_{mn} are the (randomly chosen) initial phases, are positive integers, and (for convenience) *M* has an integer square root.

As with all such algorithms, the Runge–Kutta–Fehlberg algorithm basically involves adding small increments to a function—here **x**(*t*)—where the increments are given by derivatives ((d**x**/d*t*)=**v**=∇*S*=Im (∇*Ψ*/*Ψ*)) multiplied by variable step sizes (here, a time step Δ*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, then a smaller time step 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, then 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 time scales 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 16 processors (17 evaluations of between time 0 and time 4*π*) took 9.6 CPU-hours with a nine-mode wave function. A comparable calculation with 36 modes took 542 CPU-hours. The calculations we report here are for 4, 9, 16, 25, 36, 49 and 64 modes. Previous calculations were done exclusively with either four modes (Colin & Struyve 2010) or 16 modes (Valentini & Westman 2005).

## 3. Results and discussion

We begin by verifying that exponential decay is an appropriate model for the evolution of under these conditions, as was assumed in the definition of the relaxation time scale. Figure 4 shows plots of versus *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 over this range of conditions. This was not unexpected as exponential decay was previously demonstrated by Valentini & Westman (2005), 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*

Figure 5 shows plots of against 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 time scale with the standard deviation taken to be the error bar.

In the case of four-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 figure 5). Some representative values of *τ* were 980±1600 for *ε*=4 and 100±110 for *ε*=32. This large spread in the time scale 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 (Valentini & Westman 2005), and so the initial positions of the nodes are likely to affect the time scale (these initial positions are moved around when modifying the set of initial phases *θ*_{mn}). The important point is that in larger 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, perhaps containing only 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 in Towler *et al.* (2010). 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 figure 6, therefore, we chose to plot the same data as in figure 5 but with the points for *M*=4 omitted, in an attempt to show the relationship between *τ* and *M* in the regime where the distribution of nodes may be considered roughly fixed. Can we justify this? Arguably, points with such a large error will have little effect on the curve fitting, as points are weighted according to the inverse of their error. However, a better reason for neglecting the points for *M*=4 is that the theoretical predictions are actually in terms of Δ*E*, the uncertainty in the energy of the wave function, rather than the number of modes, *M*, in the wave function. 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 summarized in table 1. Clearly, these results do *not* support the theoretical prediction (described as a ‘crude estimate’ in Valentini & Westman (2005)) 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 equation (1.9) are invalid.

The relevant arguments used to obtain the apparently incorrect estimate of the scaling are set out in Valentini (2001). They begin by defining a relaxation time scale *τ* in terms of the rate of decrease of near *t*=0 via the following formula:
3.1The second derivative is used rather than the first derivative since (as may be shown analytically (see Valentini 1992)). The -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 begins soon after *t*=0, and that in the limit the decay is not exponential. The time scale of equation (3.1) then applies only to this very short period immediately after *t*=0. In other words, although equation (3.1) may be a good estimate for the relaxation time scale in the limit , 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 equation (1.9) is the requirement for a 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 , then applying dimensional analysis to find the other dependencies. The inverse dependence of *τ* (as defined by equation (3.1)) 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, figure 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 equation (1.9), the domain of validity of which 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 time scales *τ* such that the said upper bound becomes of the order of the width *L* of the potential well. As we shall now show, the time scale will then be
3.2where is the mean energy, and .

There now follows a derivation of the upper bound on the mean displacement *δq*_{r}. The derivation is based on Valentini (2008), but with some differences that are highlighted below.

First, note that the final displacement *δq*_{r}(*t*_{f}) has modulus
3.3where is the component of the de Broglie–Bohm velocity in the *q*_{r} direction.

The equilibrium mean 〈|*δq*_{r}(*t*_{f})|〉_{eq} then satisfies
3.4The equilibrium mean speed, , is
3.5Using the fact that, for any *x*, , we have
3.6Now note that, starting from the guidance equation,
3.7where is the momentum operator conjugate to *q*_{r}, and denotes the usual quantum expectation value for the operator . The last equality follows from the relations
3.8and
3.9Since (∂|*Ψ*|/∂*q*_{r})^{2}≥0, we then have
3.10and so
3.11We also have
3.12where *W*_{r} denotes the *x*- or *y*-part of the mean Hamiltonian, with *W*_{r}∝*M*^{2}. Hence,
3.13Since *W*_{r} is time independent, and setting *t*_{i}=0 and *t*_{f}=*t*, we have
3.14Setting the right-hand side to be of order *L*, and noting that , this indeed yields the relaxation time in equation (3.2), the inverse scaling with *M* of which is in agreement with the numerical results presented above.

As mentioned, this derivation is based on that in Valentini (2008) but differs in some respects. The purpose of the analysis in Valentini (2008) 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*_{r}—for 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.

### (b) Relaxation time as a function of coarse-graining length

In figure 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 2. It is evident from these data that, although the relaxation time evidently decreases with increasing *ε*, the data do not support equation (1.9) nor are they particularly suggestive of a constant power law. A weak dependence of order is observed. Equation (1.9) predicts a power-law relation *τ*∝*ε*^{−1} while equation (3.2) predicts no dependence.

As was the case with the *M*-dependence of *τ*, the dependence on *ε* is not in fact expected to take the form in equation (1.9), for two reasons. First, the derivation of equation (1.9) is based on a definition of ‘relaxation time’ that applies only when ; it does not apply to the exponential tail of , 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 equation (1.9). Indeed the graphs in figure 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 (which must hold as we have noted), since at some point to the left of the data in any of the graphs 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 equation (3.2) is not surprising since the analysis in §3*a* and Valentini (2008) uses a different definition of the time scale. This definition does not consider the -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.

## 4. 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, the Born rule does *not* have to be introduced as a postulate of non-relativistic quantum mechanics. What is the price paid for this? We must suppose merely that particles have well-defined positions (and hence trajectories) *continuously* rather than only when a position measurement is performed.

The main technical result of this work is the emergence of the relationship *τ*∝*M*^{−1} showing the dependence of the relaxation time on the number of modes in a superposition. This result seems fairly robust under the conditions studied here. In general terms, the faster relaxation times for larger *M* are due to the greater number of free-moving nodes in the pilot wave which act as a source of vorticity and increase the chaotic nature of trajectories. Our numerical result for the scaling conflicts with the previous theoretical prediction, equation (1.9), but agrees with an alternative theoretical analysis presented here in §1*b*. As we have discussed, the assumptions made in deriving equation (1.9) were probably too restrictive for it to be of practical use. In particular, the defined time scale is relevant only close to *t*=0, and does not apply to the exponential tail of the 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 . This also differs from equation (1.9) and this is probably simply because the coarse-graining cells are too large for the derivation of equation (1.9) to be valid. It is possible that by decreasing *ε* a behaviour conforming better to equation (1.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 approximately 1 and 1000. Using natural units 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 Valentini (2008, 2010) and works cited therein) that relic non-equilibrium from the early Universe could be observed today, either directly or by its imprint on the CMB. 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 Valentini (2008, 2010). 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 (Valentini 1991*a*,*b*), 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 subensembles that obey the Born rule. The original view expressed in Valentini (1991*a*,*b*) 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 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 subensemble.)

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 CMB. 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 Valentini (2010).) 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 3*N* dimensions.

## Footnotes

↵1 In the high-energy domain, pilot-wave theory for bosons usually takes the form of a field theory, while for fermions the best model invokes the Dirac sea. These models have a fundamental preferred rest frame, with an effective Lorentz invariance emerging only in equilibrium. For recent progress, see in particular Colin (2003), Colin & Struyve (2007) and Struyve (2010).

↵2 It should be noted, however, that, in the code used by Valentini and Westman, the signs of the initial phases in the wave function were chosen with the opposite convention to that given in their text (where the latter agrees with equation (2.3)). Furthermore, in their plots of density functions, the labels on the

*x*- and*y*-axes were inadvertently exchanged.↵3 Assuming

*ρ*and |*Ψ*|^{2}are normalized, then . Since for any value of*ρ*and |*Ψ*|^{2}, with equality only when*ρ*=|*Ψ*|^{2}, then it is clear that . Thus,*H*is bounded from below by zero.

- Received October 3, 2011.
- Accepted November 1, 2011.

- This journal is © 2011 The Royal Society

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.