## Abstract

Numerical simulations indicate that the Born rule does not need to be postulated in the de Broglie–Bohm pilot-wave theory, but tends to arise dynamically (relaxation to quantum equilibrium). These simulations were done for a particle in a two-dimensional box whose wave function obeys the non-relativistic Schrödinger equation and is therefore scalar. The chaotic nature of the de Broglie–Bohm trajectories, thanks to the nodes of the wave function which yield to vortices, is crucial for a fast relaxation to quantum equilibrium. For spinors, we typically do not expect any node. However, in the case of the Dirac equation, the de Broglie–Bohm velocity field has vorticity even in the absence of nodes. This observation raises the question of the origin of relaxation to quantum equilibrium for fermions. In this article, we provide numerical evidence to show that Dirac particles also undergo relaxation, by simulating the evolution of various non-equilibrium distributions for two-dimensional systems (the two-dimensional Dirac oscillator and the Dirac particle in a two-dimensional spherical step potential).

## 1. Introduction

In the de Broglie–Bohm pilot-wave theory (de Broglie 1928; Bohm 1952*a*,*b*), each element of an ensemble (consisting of a single particle) is described by a wave function *ψ*(*t*,** x**) ( in the polar representation) but also by a position

**(**

*x**t*). The wave function always evolves according to the Schrödinger equation 1.1whereas the evolution of

**(**

*x**t*) is determined by the standard guidance equation 1.2The de Broglie–Bohm theory reproduces the predictions of standard quantum mechanics provided that the particle positions are distributed according to 1.3over the ensemble. de Broglie's dynamics (defined by equation (1.1) and equation (1.2), by contrast to Bohm's dynamics (Colin

*et al.*in preparation), which is a second-order formulation of the dynamics based on the quantum potential) ensures that this condition will hold at all times, if it holds for some initial time

*t*

_{0}. This property is referred to as equivariance (Dürr

*et al.*1992), the distribution given in equation (1.3) being the quantum equilibrium distribution (Valentini 1991

*a*,

*b*). There are (at least) four different attitudes towards this condition (and some of them are not mutually exclusive): Valentini (1991

*a*,

*b*, 1992) argues that it arises dynamically, Dürr

*et al.*(1992) argue that we happen to live in a typical universe, Holland (1993) and Bricmont (2001) take it as a third postulate, while Wiseman suggests that it can be derived as a Bayesian prior from a principle of indifference (Wiseman 2007). Quantum non-equilibrium, which leads to new physics and testable predictions of the pilot-wave theory, has been explored by Valentini (see, for instance, Valentini (2010) or Colin & Wiseman (2011) for related work).

If the possibility of quantum non-equilibrium is considered seriously, one has to explain why we do not see quantum non-equilibrium today. This is done by invoking the idea of relaxation to quantum equilibrium: if some non-equilibrium distributions existed in the past, they are quickly driven dynamically to quantum equilibrium. Various numerical simulations for two-dimensional systems (Valentini & Westman 2005; Colin & Struyve 2010; Towler *et al.* 2011) have shown that non-equilibrium distributions rapidly relax to quantum equilibrium on a coarse-grained level, provided that the wave function has enough complexity. By complexity, we mean that two trajectories originating from neighbourhood points should quickly diverge (which is referred to as chaos in a loose sense). The existence of nodes has been first connected to chaos in Frisk (1997). So when we say that the wave function needs to be complex enough, we mean that one needs to superpose a few energy modes in order to get a few nodes. Nodes are also associated to vorticity: the circulation of the velocity field around a closed curve can only be non-zero if there is a node inside the curve. More recently, a better understanding of chaos, through ‘nodal-point–X-point complex’ (where the X-point is an unstable point associated to the node), has also been gained (Efthymiopoulos *et al.* 2009).

Although the previous relaxation simulations have some bearing to non-relativistic fermions, they strictly apply to non-relativistic scalar particles. Also no simulation has ever been done for relativistic fermions. And there is a good reason to do simulations for fermions, because we typically do not expect any node for spinors, no matter how many energy modes are being superposed. The reason is the following. For scalar particles, we have nodes where and . In two dimensions, nodes are located at the intersection of the curves defined by the two previous equations and we typically have a few nodal points (in three dimensions, we would have nodal lines). For a positive-energy spin- fermion, described by a Pauli spinor , nodes are located at the intersection of four curves (or surfaces): , , and . Therefore, we typically do not have any node in two or three dimensions. This raises the question of relaxation to quantum equilibrium. At first, it seems problematic for spinors.

But the velocity field for spinors (for Dirac particles for instance) is essentially different from the standard velocity-field for non-relativistic scalar particles: it has vorticity even in the absence of nodes. If vorticity is the crucial ingredient for relaxation to quantum equilibrium, we expect relaxation to quantum equilibrium for Dirac particles too. Actually, a previous work by Colin & Struyve (2010) gives some support to that idea. There the evolution of a non-equilibrium distribution in a two-dimensional square box is simulated, the wave function being nodeless (except at the boundaries). Relaxation is very poor in that case, except if ‘artificial vortices’ are added, such that the alternative (and non-standard) velocity field exhibits vorticity even in the absence of nodes. This is always possible and is referred to as the non-uniqueness problem of the de Broglie–Bohm theory (Deotto & Ghirardi 1998; Wiseman 2007). The example of the nodeless wave function with artificial vortices can be a toy-model for spinors; however, one can argue with reason that these non-standard velocity fields are unnatural and that this example lacks features which are inherent to spinors (for instance, the spinorial density cannot vanish at the boundaries).

That is the main question we will address in this paper. In order to do that, we will consider 2 two-dimensional systems for ‘confined’ Dirac particles: the first one is the so-called Dirac oscillator, the second one is a Dirac particle confined by a spherical step potential. In each case, we study a superposition of positive-energy solutions. We will show that the first system, in which eight modes are superposed, exhibit fast relaxation, while the second system, in which six modes are superposed, also undergoes relaxation to quantum equilibrium (although it is a slower one). In the second system, we also give an example of a non-trivial spinor for which there will never be complete relaxation.

This article is organized in the following way. Section 2 contains the basic details about the Dirac equation in 3+1 and 2+1 dimensional space–times, as well as the corresponding pilot-wave theory. In §§3 and 4, we study relaxation to quantum equilibrium for Dirac spinors. We conclude in §5.

We use units in which .

## 2. Pilot-wave theory for a Dirac fermion

In this section, we recall the essential properties of the Dirac equation in space–times of dimension 4 and 3, and the corresponding pilot-wave theory.

### (a) 3+1-dimensional space–time

The Dirac equation is
2.1where the matrices *α*_{j} and *β* must be Hermitian and must satisfy the relations {*α*_{j},*α*_{k}}=2*δ*_{jk}, *β*^{2}=1 and {*α*_{j},*β*}=0.

The Dirac equation can be rewritten in a covariant form by introducing the matrices *γ*^{0}=*β* and *γ*^{j}=*βα*_{j}:
2.2The *γ*-matrices satisfy the relations {*γ*^{μ},*γ*^{ν}}=2*η*^{μν}, where *η*^{μν}=diag(1,−1,−1,−1) (Clifford algebra). The Pauli–Dirac representation
2.3and the Weyl (or chiral) representation
2.4are two common choices of representation for the *γ*-matrices. The conserved four current is given by
2.5where .

The positive and negative-energy plane-wave solutions are denoted by *u*(** p**)e

^{−iE(p)t+ip⋅x}and

*v*(

**)e**

*p*^{iE(p)t+ip⋅x}, where . In the Weyl representation, if one introduces the right-handed and left-handed eigenstates of helicity

*χ*

_{R}(

**) and**

*p**χ*

_{L}(

**), the plane-wave solutions can be given by 2.6and 2.7**

*p*In the corresponding pilot-wave theory (first proposed by Bohm in a reply to Takabayasi (Bohm 1953)), the particle is described by a four spinor *ψ*(*t*,** x**) together with its position

**(**

*x**t*). The position of the Dirac particle evolves according to the guidance equation 2.8Again the choice of the guidance equation ensures that the relation 2.9will hold for any time

*t*provided that it holds for some initial time.

Any positive-energy plane-wave solution of momentum ** p**, whether right-handed or left-handed, moves with velocity

**/**

*p**E*, whereas negative-energy plane-wave solutions move in the opposite direction of momentum.

Figure 1*a* shows the kind of trajectories predicted by the theory. The guiding Dirac spinor in this case is a superposition of three positive-energy right-handed plane-wave solutions
2.10where
2.11and *u*_{R} is defined at equation (2.6).

### (b) 2+1-dimensional space–time

In a 2+1-dimensional space–time, the Dirac equation is
2.12We can take *α*_{1}=*σ*_{1}, *α*_{2}=*σ*_{2} and *β*=*σ*_{3}.

The covariant form is
2.13where the *γ*-matrices are defined by
2.14

For the plane-wave solutions, we have 2.15and 2.16which are, respectively, positive and negative-energy solutions.

In figure 1*b*, we plot a few trajectories of particles guided by the spinor
2.17where *q*_{1}=(1,0), *q*_{2}=(−1,−2) and *q*_{3}=(1,−1), for different values of the mass.

The circulation of the velocity field (where ** v**=(

*ψ*

^{†}

*σ*

_{1}

*ψ*,

*ψ*

^{†}

*σ*

_{2}

*ψ*)/

*ψ*

^{†}

*ψ*) can be non-zero, indicating vorticity, even in the absence of nodes. It results from an application of Stokes' theorem and from the fact that the curl of the Dirac velocity field can be non-zero (on the other hand, the curl of the standard non-relativistic velocity field is always zero except at nodes where the velocity field is not defined).

## 3. Relaxation to quantum equilibrium for the two-dimensional Dirac oscillator

The equation for the Dirac oscillator is obtained by substituting ** p** by

**−**

*p**imωβ*

**in the free Dirac Hamiltonian (Moshinsky & Szczepaniak 1989): 3.1This system reduces to a standard oscillator plus a strong spin-orbit coupling term in the non-relativistic limit and some of its generalizations (Castaños**

*r**et al.*1991; Kukulin

*et al.*1991) are relevant to the study of quarks confinement.

The two-dimensional case 3.2with Hamiltonian 3.3was treated by Villalba (1994). We have verified Villalba's construction and we arrive at the same energy eigenstates expect that a few intermediate equations are different, which is why we include the derivation here. Also, it will be useful for the case of the Dirac particle confined by the step potential.

### (a) Energy eigenstates

If the Cartesian coordinates are replaced by polar coordinates, i*H* becomes:
3.4The following substitution
3.5is then made in equation (3.2). We first evaluate the action of the operator ** α**⋅

**(first matrix in equation (3.4)) on the spinor: 3.6and 3.7Thanks to equations (3.4)–(3.7), equation (3.2) becomes 3.8and 3.9Now we look for energy eigenstates: we make the replacement 3.10where**

*p**k*is half-integer (and equal to the total angular momentum), in order to find 3.11and 3.12from which the following equation for can be deduced: 3.13Now we perform the following change of variables

*ξ*=

*mωr*

^{2}, which implies that 3.14Equation (3.13) becomes 3.15where Δ

^{2}=(

*E*

^{2}−

*m*

^{2})/

*mω*. If we make the following replacement 3.16in equation (3.15), we obtain the equation 3.17It reduces to an equation for generalized Laguerre polynomials 3.18provided that

*α*=

*k*/2 (or

*α*=1/2−

*k*/2),

*μ*=

*k*−1/2 (or −

*k*+1/2) and 3.19The choice of

*α*(which in turn fixes

*μ*and Δ

^{2}) is fixed by the sign of

*k*, because the wave function needs to be regular at the origin (for example, if

*k*is positive, equation (3.16) forces us to choose

*α*=

*k*/2). can then be obtained from equation (3.12).

We can always choose units in which , *c*=1 and *m*=1. The frequency *ω* remains a free parameter. The following table contains the eigenstates of lowest energy for the case *ω*=1.
3.20and
3.21

### (b) Simulations

For the wave function, we consider a superposition of the eight previous energy eigenstates with equal weights and phases given in the following table:
3.22We consider five different non-equilibrium distributions at time *t*=0. The first one, denoted by *ρ*_{0}(*t*=0,*r*,*θ*), is defined by
3.23if *r*≤*R*_{0}, and zero otherwise. The four remaining distributions (*ρ*_{1},*ρ*_{2},*ρ*_{3} and *ρ*_{4}) are built from *ρ*_{0}, except that they are, respectively, centred at *x*_{1}=(2,0), *x*_{2}=(0,2), *x*_{3}=(−2,0) and *x*_{4}=(0,−2) and that *R*_{0} is substituted by *R*=2. Explicitly, we have that:
3.24if , and zero, otherwise.

The evolution of the non-equilibrium distribution can be computed using a method proposed by Valentini & Westman (2005) and also used later in Colin & Struyve (2010) and Towler *et al.* (2011). This method relies on the fact that the following ratio
3.25is preserved along a trajectory. Therefore, if we want to compute *ρ*(*t*,** x**), we first ‘backtrack’ the initial position (that is, we find the initial condition (

*t*

_{0},

*x*_{0}) which gives (

*t*,

**) as final position). Then, we use equation (3.25) in order to find that 3.26For this simulation, we have used a lattice of 2048×2048 points, uniformly distributed inside the box [−5,5]×[−5,5], of coordinates 3.27with**

*x**k*,

*l*∈{0,1,2,…,2047}. Note that all the lattice points are inside the box (none of them lies on the boundary of the box), and that is why 10/4098 appears in the expression (and not 10/2048). Each lattice point (representing the final position of a particle) is backtracked to its initial position by solving the differential equation 3.28numerically using the Runge–Kutta–Fehlberg (RKF) method (Press

*et al.*1992). We require a precision of 0.001 on the backtracked positions (this means that the distance between two backtracked positions resulting from two successive iterations of the RKF method with absolute tolerance parameters

*ϵ*and

*ϵ*/10 should be smaller than 0.001). Lattice points that give an insufficient precision for the chosen minimal value of

*ϵ*, or cannot be backtracked before a given maximum number of iterations (for instance 100 000), are discarded. Hereafter, we refer to these lattice points as bad ones (or good ones if the converse is true).

The next step is the coarse-graining, which is done by dividing the box [−5,5]×[−5,5] in a certain number of non-overlapping coarse-graining (CG) cells, and averaging the densities corresponding to the lattice points contained in each CG cell. If we divide the box in 32×32 non-overlapping CG cells, the mean percentage of good lattice points per CG cell is 98.47 per cent and the worst cell has 85.06 per cent.

For the plots, we have actually used a smooth coarse-graining of the densities (the same that was used in Colin & Struyve (2010)) which is similar to a standard coarse-graining except that it is defined over overlapping CG cells. In this case, all the overlapping CG cells can be obtained from a square cell of side 10/16 located in the lower left corner by translating it by in the up-direction, and by in the right-direction, where *m*,*n*∈{0,120}. The smooth coarse-graining is denoted by
3.29where ** x** is in fact the centre of an overlapping CG cell.

The evolutions of the non-equilibrium distributions *ρ*_{0} and *ρ*_{1} are illustrated by density plots and surface plots in figure 2.

### (c) Discussion

For that superposition of eight modes with positive energy, the relaxation to quantum equilibrium is rather fast. In these units, it takes about Δ*t*=100 for the non-equilibrium distribution to resemble |*ψ*|^{2}, which is about one hundred times the Compton wavelength of the particle divided by the speed of light. In standard units, it is about 10^{−18} s if we take the Compton wavelength of the electron, which is of the order of 10^{−12} m.

## 4. Relaxation to quantum equilibrium for a Dirac particle in a two-dimensional spherical step potential

We consider a Dirac particle confined by a spherical step potential: 4.1The corresponding three-dimensional case is analysed in Greiner (1990).

### (a) Energy eigenstates

It is similar to the previous case: we can start from equation (3.13), replace *E* by *E*+|*V* | (where |*V* | is equal to |*V* _{0}| or equal to zero, depending on whether we consider the interior solution or the exterior one) and set *ω*=0. Then we have that
4.2where *k* is half-integer (we define ).

We distinguish two cases:

*Case 1*. *κ*^{2}=(*E*+|*V* |)^{2}−*m*^{2}>0. In that case, the equation to solve is
4.3If we introduce *s*=*rκ* and , we find the following equation:
4.4which admits Bessel functions of the first and second kind as solutions (*J*_{k−1/2}(*s*) and *Y* _{k−1/2}(*s*)). If we put all the pieces together, we have that
4.5and
4.6For , we have that
For *ψ*_{2},
4.7

*Case 2*. −*κ*^{2}=(*E*+|*V* |)^{2}−*m*^{2}<0. In that case, the equation is
4.8If we introduce again *s*=*rκ* and , we find the following equation:
4.9It admits modified Bessel functions of the first and second kind as solutions (denoted by *I*_{k−1/2} and *K*_{k−1/2}).

We find that 4.10and that 4.11

The next step is to choose some numerical values for the potential *V* _{0} and the radius of the spherical step potential and find the energy eigenvalues numerically.

### (b) Simulations

We choose *m*=1, *R*′=5 and |*V* _{0}|=1. We are going to superpose energy eigenstates whose eigenvalues *E* satisfy the relations
4.12Therefore, we take

— equations (4.6) and (4.7) for the interior solution (and we set

*β*=0 for the solutions to be regular at the origin),— equations (4.10) and (4.11) for the exterior solution (and we set

*α*′=0 for the solutions to vanish at ).

The next step is to find the energy eigenvalues numerically, for different values of *k*, by imposing that
4.13and by looking in the domain defined by equation (4.12). The last step is to find *β*′ by requiring that
4.14Overall, these two conditions amount to the following boundary condition: the interior solution for the spinor should coincide with the exterior one at *r*=*R*′.

The guiding spinor used for the simulation is a superposition of six positive-energy eigenstates with different values of *k*. The details about the spinor can be found in appendix A.

For the relaxation simulations, we have considered the evolution of the five non-equilibrium distributions for *t*=50 and *t*=100. The lattice consists of 1024×1024 points distributed uniformly in the box [−3,3]×[3,3] and the precision on backtracked trajectories is 0.005. If the box is divided into 32×32 CG cells, the mean percentage of good trajectories at *t*=100 is 93.22 per cent and the worst cell has a percentage of 67.58 per cent. Some of these results are illustrated in figure 3*a*,*b*.

### (c) Discussion

While the relaxation to quantum equilibrium is slower in this case (compared with the Dirac oscillator case), it must be pointed out that the CG length is smaller here (3/5 of the CG length used in the previous case). Also, if the relaxation time depends on the energy spreading, we would expect a longer relaxation time in this case anyway. Finally, the Dirac oscillator has a strong spin–orbit coupling (Kukulin *et al.* 1991) which may speed up the relaxation.

Another interesting example is the case of a wave function with only positive values of *k* (which means that the direction of rotation does not change). In that case, non-equilibrium distributions will not relax (figure 4). In order to understand this, we plot a few backtracked trajectories (last plot of figure 4) for a superposition of the three modes with positive values of *k*. As we can see, trajectories from the inner core (which is a small ball centred at the origin) are backtracked to the inner core, which is a bad thing for relaxation. One might argue that only the modes with *J*_{0}(*x*) are dominant in that region (these are the eigenstates with *k*=1/2 and *k*=3/2). We have only tested that the inner core is backtracked to the inner core for a superposition of three modes with positive values of *k* but it might be a general feature of many superpositions with positive values of *k*. For instance, if we simply add one mode with *k*=−1/2, the inner core is not disconnected from the outer core anymore: trajectories backtracked from the inner core can end up in the outer core.

## 5. Conclusion

We have considered 2 two-dimensional systems (the Dirac oscillator and the Dirac particle in a spherical step potential) and we have done numerical simulations which show an efficient relaxation to quantum equilibrium for Dirac fermions, provided that enough modes are superposed in order to get sufficient complexity in the dynamics. The relaxation takes place despite the absence of nodes and it seems natural to attribute it to the intrinsic vorticity of the de Broglie–Bohm Dirac velocity field (although the roles played by X-points and quasi-nodes should be further studied). For the second system, we have also given an example of a non-trivial spinor (with positive values of *k*) for which there will never be complete relaxation.

This last example does not ruin the idea of relaxation. Firstly, it is a two-dimensional system. But secondly (and more importantly), even in a two-dimensional universe, this system does not sit there on its own: it starts to interact with other systems, resulting in a more complex system, where relaxation can take place. However, this example is interesting with respect to relic non-equilibrium particles (Valentini 2007), which are hypothetical non-equilibrium particles from the very early universe that would still be in non-equilibrium today. The possibility of finding such systems arises because the spatial expansion in the very early universe stretches non-equilibrium length-scales and competes with the natural process of relaxation. Then, if these non-equilibrium particles decouple very early and if they are still in non-equilibrium at the time of decoupling, they cannot be driven to quantum equilibrium by interaction with other systems. Furthermore, if they are in an energy eigenstate, the natural relaxation does not take place and the non-equilibrium can be preserved up to the present day. The last example shows that non-equilibrium could also be preserved when the spinor is more complex than an energy eigenstate.

As possible future work, one can study

— the existence and roles of X-points and quasi-nodes for the Dirac velocity field,

— the influence of spatial expansion (which is relevant for the early universe and results in a stretching of the non-equilibrium length-scales),

— the role of mass, since it is not fundamental and beables should really be attributed to massless fermions (discussed in Colin & Wiseman (2011)),

— relaxation timescales, via the H-function (like in Towler

*et al.*(2011)) and— the generalization to 3+1-dimensional systems, in particular spinning systems for which the rotation direction does not change.

Also the Dirac equation is not the only relativistic equation. For instance, it is still an open question whether neutrinos are Dirac or Majorana fermions (a Majorana particle being its own anti-particle). The Majorana equation is not only relevant for particle physics but also to quantum information and condensed matter research (an overview of the wide relevance of the Majorana equation is given in Wilczek (2009)). From the point of view of the pilot-wave theory, the Majorana equation is also very particular. Indeed, in a forthcoming work, we will show that the Majorana equation predicts luminal motion for the beable, although the Majorana solutions involve the mass. In the same paper, we study relaxation to quantum equilibrium for systems governed by the Majorana equation.

## Acknowledgements

I thank Howard Wiseman and a referee for comments and suggestions on a previous manuscript. I acknowledge financial support from a Perimeter Institute Australia Foundations post-doctoral research fellowship. This work was also supported by the Australian Research Council Discovery Project DP0880657, ‘Decoherence, Time-Asymmetry and the Bohmian View on the Quantum World’. Research at Perimeter Institute is funded by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. The relaxation simulation for the Dirac oscillator was performed on the *Griffith University HPC Cluster, V20z*.

## Appendix A. Guiding spinors for the second system

We take |*V* _{0}|=1, *m*=1 and *R*′=5.

For an energy eigenstate, the interior solution is
A1and
A2where , while the exterior solution (*r*≥*R*′) is
A3and
A4with .

The energy eigenvalues and the parameters *β*′ are found by matching *ψ*_{1} and *ψ*_{2} for the interior and exterior solutions:
A5Finally, we take random phases e^{iϕ}:
A6The superposition of six modes is obtained by summing the corresponding eigenstates multiplied by the phase coefficients. In the end, the spinor is normalized numerically. The superposition of three (respectively, four) modes is obtained by superposing only the first three (respectively, four) eigenstates.

- Received September 11, 2011.
- Accepted November 14, 2011.

- This journal is © 2011 The Royal Society