Master equations for Wigner functions with spontaneous collapse and their relation to thermodynamic irreversibility

Abstract

Wigner functions, allowing for a reformulation of quantum mechanics in phase space, are of central importance for the study of the quantum-classical transition. A full understanding of the quantum-classical transition, however, also requires an explanation for the absence of macroscopic superpositions to solve the quantum measurement problem. Stochastic reformulations of quantum mechanics based on spontaneous collapses of the wavefunction are a popular approach to this issue. In this article, we derive the dynamic equations for the four most important spontaneous collapse models—Ghirardi–Rimini–Weber (GRW) theory, continuous spontaneous localization (CSL) model, Diósi-Penrose model, and dissipative GRW model—in the Wigner framework. The resulting master equations are approximated by Fokker–Planck equations. Moreover, we use the phase-space form of GRW theory to test, via molecular dynamics simulations, David Albert’s suggestion that the stochasticity induced by spontaneous collapses is responsible for the emergence of thermodynamic irreversibility. The simulations show that, for initial conditions leading to anti-thermodynamic behavior in the classical case, GRW-type perturbations do not lead to thermodynamic behavior. Consequently, the GRW-based equilibration mechanism proposed by Albert is not observed.

1 Introduction

Wigner functions provide a description of quantum systems in phase space. Since their initial development [1], they have found a significant number of applications [2,3,4,5] in fields such as atomic physics [6], quantum optics [7,8,9], visualization of quantum effects [10, 11], computational electronics [12,13,14], and solid-state theory [15,16,17,18,19]. Besides these practical aspects, they are also of interest for fundamental questions in quantum mechanics such as the theory of quantum chaos [20]. Moreover, they have a natural connection to statistical mechanics [21, 22], including nonequilibrium descriptions such as the Boltzmann equation [23,24,25] and the Fokker-Planck equation [26,27,28]. In particular, it has been suggested that Wigner functions might be useful in understanding the origin of thermodynamic irreversibility [29].

Since they allow to describe quantum mechanics in a way that has strong formal analogies to classical mechanics (Wigner functions are defined on phase space and governed by a dynamic equation that reduces to the classical Liouville equation for the phase-space density in the limit $\hbar \rightarrow 0$), an important application of Wigner functions is the study of the quantum-classical transition [30, 31]. Wigner functions allow to develop a connection between classical and quantum models for liquid crystals [32] and to study mixed quantum-classical dynamics [33,34,35]. However, recovering the Liouville equation is not sufficient for recovering classical mechanics. Quantum systems can be in superpositions, whereas macroscopic superpositions (such as a cat being dead and alive at the same time [36]) are never observed. This is commonly attributed to a “collapse of the wavefunction” that takes place when a measurement is performed. However, this approach is often seen as unsatisfactory given that such a collapse is not allowed by the Schrödinger equation, meaning that one has to make the strange assumption that the Schrödinger equation does not apply to measurements. This issue is known as the measurement problem [37].

A popular approach to the measurement problem is to propose a stochastic modification of quantum mechanics in which collapses occur spontaneously [38,39,40,41]. Spontaneous collapse models make empirical predictions that differ slightly from those of “standard” quantum mechanics and can thus be tested experimentally [42, 43]. Such experiments have become an important field of research in recent years [44,45,46]. Incorporating such spontaneous collapses in the theory of Wigner functions is desirable for at least four reasons: First, it allows to solve the problem of the classical-quantum transition in the Wigner framework also regarding the measurement problem. Second, if collapse models should be confirmed experimentally, the Wigner formalism should of course allow to describe them in order to remain a powerful alternative to Hilbert-space quantum mechanics. Third, Wigner functions are very useful for the visualization of collapse effects, as has been exploited for both spontaneous [47, 48] and environment-induced [49] collapses. Fourth, modern quantum technology—a field where Wigner functions are of high importance [2]—may be relevant for experimental tests of spontaneous collapse models. For example, it has been discussed whether SQUIDs (superconducting quantum interference devices) can be used for such a test [50, 51].

In addition, spontaneous collapse models might provide a solution to another fundamental problem of physics, namely to the aforementioned problem of thermodynamic irreversibility [52, 53]. This problem (which actually includes a variety of subproblems [53]) is concerned with the fact that macroscopic thermodynamics has a clear arrow of time associated with the monotonous increase of entropy, whereas the microscopic laws of (classical or quantum) mechanics are invariant under time reversal. David Albert [54,55,56] has suggested that spontaneous collapses of the wavefunction might be responsible for the irreversible approach to equilibrium. Although it has received significant attention in philosophy of physics [57,58,59,60,61,62], rigorous physical tests of this proposal are lacking at present.

In this work, we provide a reformulation of the four most important spontaneous collapse models, the Ghirardi-Rimini-Weber (GRW) theory [38], the continuous spontaneous localization (CSL) model [39, 63], the Diósi-Penrose model [64,65,66,67], and the dissipative GRW model [68, 69], based on Wigner functions. We then provide an analytical and numerical analysis of Albert’s proposal, thereby investigating whether GRW theory might actually contribute to the solution of the problem of irreversibility. The mechanism suggested by Albert (deviations resulting from stochastic perturbations lead from anti-thermodynamic to thermodynamic behavior) is not observed in the simulations.

2 Structure of this article

Since this article is rather long and involves results of interest to a variety of audiences (ranging from quantum physicists to philosophers of physics), we here briefly summarize its structure to enable readers interested only in a particular topic to quickly find the relevant section.

The article consists of two main parts. Although both contain important results, a reader only interested in one of them may skip the other one. The first part deals with the foundations of quantum mechanics. In Sect. 3, we provide a brief introduction to Wigner functions. The following four sections each explain a particular collapse model and derive the equation of motion that it implies for the Wigner function. This is done in Sect. 4 for the GRW theory [with the main results being Eqs. (13) and (15)], in Sect. 5 for the CSL model [with the main result being Eq. (28)], in Sect. 6 for the Diósi-Penrose model [with the main result being Eq. (39)], and in Sect. 7 for the dissipative GRW model [with the main results being Eqs. (42) and (44)].

The second part deals with the foundations of statistical mechanics, in particular with Albert’s suggestion that GRW theory might explain thermodynamic irreversibility. This idea is introduced and discussed in Sect. 8. In Sect. 9, decoherence-based alternatives are presented. A numerical test of the GRW-based explanation of irreversibility is provided in Sect. 10. The main result of the second part are the simulations shown in Fig. 1, which disagree with Albert’s proposal. We conclude in Sect. 11.

3 Wigner functions

In “standard” quantum mechanics, a system is described by the density operator¹ (also known as “density matrix” or “statistical operator”) ${\hat{\rho }}$ whose time evolution is given by the Liouville-von-Neumann equation

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}{\hat{\rho }}= -\frac{\mathrm {i}}{\hbar }[{\hat{H}},{\hat{\rho }}] \end{aligned}$$

(1)

with time t, imaginary unit $\mathrm {i}$, reduced Planck constant $\hbar$, commutator $[\cdot ,\cdot ]$, and Hamiltonian ${\hat{H}}$. Furthermore, an observable is represented by a self-adjoint operator ${\hat{A}}$.

An alternative to the density-operator-based formalism is to use Wigner functions [1]. Here, the state of the system is described by a function W(x, p) (Wigner function) that depends on the phase-space coordinates x and p. (In this work, we only require Wigner functions depending on position x and momentum p. However, the approach is much more general and can also be applied to systems with other degrees of freedom, such as spin [32, 70].) The Wigner function is a generalization of a classical phase-space distribution function, from which it differs in that it can also take negative values [71]. For a one-particle system in one dimension, the Wigner function is defined as [72]

$$\begin{aligned} W(x,p) = \frac{1}{2\pi \hbar }\int \hbox {d}{x'}\bigg \langle {x-\frac{1}{2}x' \bigg |{\hat{\rho }}\bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}. \end{aligned}$$

(2)

Throughout this work, the dependence of the Wigner function W and density operator ${\hat{\rho }}$ on time t is not written explicitly to improve readability. For N particles in three dimensions, Eq. (2) generalizes to [73]

$$\begin{aligned}&W(\vec {r}_1,\dotsc ,\vec {r}_N,\vec {p}_1,\dotsc ,\vec {p}_N) \nonumber \\&= \frac{1}{(2\pi \hbar )^{3N}}\int {\hbox {d}^3r_1'}\cdots \int {\hbox {d}^3r_N'} e^{\frac{\mathrm {i}}{\hbar }\sum _{i=1}^{N}\vec {r}_i' \cdot \vec {p}_i} \nonumber \\&\quad \, \bigg \langle {\vec {r}_1-\frac{1}{2}\vec {r}_1',\dotsc ,\vec {r}_N - \frac{1}{2}\vec {r}_N' \bigg |{\hat{\rho }}\bigg |\vec {r}_1+\frac{1}{2}\vec {r}_1',\dotsc , \vec {r}_N+\frac{1}{2}\vec {r}_N'\bigg \rangle }. \end{aligned}$$

(3)

Equations (2) and (3) can not only be applied to the density operator ${\hat{\rho }}$, but to an arbitrary operator ${\hat{A}}$. This gives its so-called Weyl symbol A(x, p).² The expectation value $\langle {\hat{A}}\rangle$ of the observable ${\hat{A}}$, which in “standard” quantum mechanics is given by

$$\begin{aligned} \langle {\hat{A}}\rangle = {\text {Tr}}({\hat{A}}{\hat{\rho }}) \end{aligned}$$

(4)

with the quantum-mechanical trace ${\text {Tr}}$, can then be calculated as [72]

$$\begin{aligned} \langle {\hat{A}}\rangle = \int {\hbox {d}x}\int {\hbox {d}p}A(x,p)W(x,p) \end{aligned}$$

(5)

just as in classical mechanics.

We now discuss the dynamics of Wigner functions. First, we introduce the star product [74]: When transforming to phase space, one has to replace the product ${\hat{A}}{\hat{B}}$ of two Hilbert space operators ${\hat{A}}$ and ${\hat{B}}$ by their star product $A \star B$, which is defined as [75]

$$\begin{aligned} A(x,p)\star B(x,p) = A(x,p)\exp \!\bigg (\frac{\mathrm {i}\hbar }{2}(\overleftarrow{\partial _x}\overrightarrow{\partial _p} -\overleftarrow{\partial _p}\overrightarrow{\partial _x})\bigg )B(x,p). \end{aligned}$$

(6)

Here, the arrows indicate in which direction the derivatives act. The star product (6) allows to introduce the Moyal bracket [74, 76]

$$\begin{aligned} \begin{aligned}&\{A(x,p),B(x,p)\}_\star \\&= -{\frac{\mathrm {i}}{\hbar }}(A(x,p) \star B(x,p) - B(x,p)\star A(x,p))\\&= \frac{2}{\hbar } A(x,p) \sin \!\bigg (\frac{ \hbar }{2}(\overleftarrow{\partial _x}\overrightarrow{\partial _p} -\overleftarrow{\partial _p}\overrightarrow{\partial _x})\bigg )B(x,p), \end{aligned} \end{aligned}$$

(7)

which is the Wigner equivalent of the commutator. When neglecting terms of order $\hbar ^2$, the Moyal bracket (7) reduces to the classical Poisson bracket. Combining Eqs. (1) and (7) and taking into account the correspondence between commutators and Moyal brackets leads to the dynamic equation for the time evolution of Wigner functions

$$\begin{aligned} \partial _t W(x,p) = \{H(x,p),W(x,p)\}_\star , \end{aligned}$$

(8)

where H is the Weyl symbol of the Hamiltonian ${\hat{H}}$. This Weyl symbol is given by $H(x,p) = p^2/(2m) + U(x)$ with the particle mass m and potential U [71]. Equation (8) reduces to the Liouville equation of classical mechanics if terms of order $\hbar$ are neglected or if the potential is at most a second-order polynomial in x [72]. (Some difficulties of the limit $\hbar \rightarrow 0$ are discussed in Ref. [71]. In particular, the classical limit cannot be taken in this simple way if the Hamiltonian generates chaotic classical dynamics [77].) An overview over theory and applications of Wigner functions is given in Refs. [2, 4, 5, 73, 78].

4 GRW theory

The formalism introduced in Sect. 3 is of interest for the Ghirardi–Rimini–Weber (GRW) theory [38], which is a modified form of quantum mechanics. It allows to deal with the famous quantum-mechanical measurement problem [37, 79,80,81,82]: If one performs a measurement on a quantum system that is in a superposition, then, assuming that both system and measurement apparatus obey the (linear, deterministic, and unitary) Schrödinger equation, the apparatus will, at the end, be in a superposition too. The quantum-mechanical time evolution inevitably leads to macroscopic superpositions, as is well known from the thought experiment on Schrödinger’s cat [36]. However, this is not what we observe. A measurement on a system in a superposition produces a definite outcome. Therefore, one typically postulates a “collapse of the wavefunction”, which is a sudden transition to an eigenstate of the observable one has measured as a consequence of measurement. This postulate, however, has no basis in the deterministic Schrödinger equation. A variety of solutions to this problem have been proposed, including hidden variables (“Bohmian mechanics” [83, 84]), superdeterminism [85, 86], modal interpretations [87], and the existence of many worlds in which all outcomes are realized (“Everett interpretation”, also known as “many-worlds interpretation” [81, 88]). It is sometimes claimed that decoherence effects (see Sect. 9) provide a solution to the measurement problem. This, however, is wrong, since decoherence alone does not determine which one of two possible measurement outcomes occurs [79, 89].

The GRW theory [38] (reviewed in Refs. [40, 42, 90]) solves the measurement problem by assuming that the collapse of the wavefunction is not a consequence of measurements, but something that happens spontaneously. For a system with N particles, the deterministic dynamics given by Eq. (1) is replaced by the more general master equation [38]

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}{\hat{\rho }}= - \frac{\mathrm {i}}{\hbar }[{\hat{H}},{\hat{\rho }}]-\sum _{i=1}^{N}\lambda _i ({\hat{\rho }}-{\hat{T}}_i({\hat{\rho }})), \end{aligned}$$

(9)

where $\lambda _i$ is the localization rate for the i-th particle. In this work, we assume that $\lambda _i = \lambda$ with a constant $\lambda$ for all particles. Since the localization rate is often assumed to depend on the particle mass [91, 92], this can be motivated by the assumption that the particles belong to the same species. ${\hat{T}}_i$ is the operator

$$\begin{aligned} {\hat{T}}_i({\hat{\rho }})= \sqrt{\frac{\alpha }{\pi }}\int _{-\infty }^{\infty }{\mathrm{d}}{x} e^{-\frac{\alpha }{2}({\hat{x}}_i - x)^2}{\hat{\rho }}\, e^{-\frac{\alpha }{2}({\hat{x}}_i - x)^2} \end{aligned}$$

(10)

with the position operator for the i-th particle ${\hat{x}}_i$ and a constant $\alpha$ determining the width of the collapsed function. It is easily shown that Eq. (9) is trace-preserving [38]. The idea of Eq. (9) is that “spontaneous collapse of the wavefunction” simply means that the particle’s wavefunction is, at a certain rate $\lambda$, stochastically multiplied by a Gaussian. The parameters are tuned in such a way that these collapses are very rare for small systems (such as quantum particles), but very common for large systems (such as a measurement apparatus). This explains why superpositions are observed in microscopic, but not in macroscopic systems.

The GRW theory (as well as other approaches to the measurement problem) is conceptually very important for a particular application of Wigner functions, namely the study of the classical limit. As discussed in Sect. 3, the governing equation (8) for Wigner functions reduces to the classical Liouville equation for $\hbar \rightarrow 0$. The measurement problem shows that recovering the Liouville equation is not sufficient to recover the observations of classical mechanics: If a cat is in a superposition of “being alive” and “being dead”, the Liouville dynamics will not change this (a superposition of a dead and a living cat evolving classically remains a superposition of a dead and a living cat). What we require is a mechanism that explains why such superpositions are not observed in the macroscopic world. The GRW theory provides such a mechanism.

A typical problem of spontaneous collapse models is that energy is not conserved. In fact, due to the noise, it increases at a constant rate. Although this effect is extremely small (the temperature of an ideal GRW gas increases by about $10^{-15}$ K per year [38]), it is still an undesirable feature. To solve this problem, dissipative extensions of the GRW theory [68], CSL model [69], and Diósi-Penrose model [93] have been derived. They ensure that the energy remains finite during the time evolution of the system. These approaches are discussed in Sect. 7. Further extensions of GRW theory (not considered in this work) are models with colored rather than white noise for the stochastic perturbations of the wavefunction [94] and relativistic extensions [95].

We now transform the governing equation (9) of the GRW theory to the Wigner formalism. We first demonstrate this derivation for the case of one particle in one spatial dimension to keep the notation compact, and then generalize to the case of N particles in three dimensions which is more relevant for statistical mechanics.

First, we use the fact that [38]

$$\begin{aligned} \langle x'|{\hat{T}}({\hat{\rho }})|x''\rangle = e^{-\frac{\alpha }{4}(x'-x'')^2}\langle x'|{\hat{\rho }}|x''\rangle , \end{aligned}$$

(11)

where $|x\rangle$ is a position eigenstate with eigenvalue x. The operator ${\hat{T}}({\hat{\rho }})$ is subject to the same transformation rule (2) as the density operator ${\hat{\rho }}$ itself. Using Eq. (11), we find

$$\begin{aligned} T(W(x,p))&=\frac{1}{2\pi \hbar }\nonumber \\&\quad \int {\hbox {d}x'}\bigg \langle {x-\frac{1}{2}x'\bigg |{\hat{T}}({\hat{\rho }})\bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}\nonumber \\&= \sqrt{\frac{\alpha }{\pi }}\int {\hbox {d}x''}\frac{1}{2\pi \hbar }\nonumber \\&\quad \int {\hbox {d}x'}\bigg \langle {x-\frac{1}{2}x'\bigg |e^{-\frac{\alpha }{2}({\hat{x}} - x'')^2}{\hat{\rho }}\, e^{-\frac{\alpha }{2}({\hat{x}} - x'')^2}\bigg |x+\frac{1}{2} x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}\nonumber \\&= \frac{1}{2\pi \hbar }\int {\hbox {d}x'}e^{-\frac{\alpha }{4}((x-\frac{1}{2}x')-(x+\frac{1}{2}x'))^2}\nonumber \\&\quad \bigg \langle {x-\frac{1}{2}x' \bigg |{\hat{\rho }} \bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}\nonumber \\&= \frac{1}{2\pi \hbar }\int {\hbox {d}x'}e^{-\frac{\alpha }{4}(x')^2}\nonumber \\&\quad \bigg \langle {x-\frac{1}{2}x'\bigg |{\hat{\rho }} \bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}\nonumber \\&= \frac{1}{2\pi \hbar }\sqrt{\frac{1}{\pi \alpha \hbar ^2}}\int {\hbox {d}x'}\int {\hbox {d}p'}e^{-\frac{1}{\alpha \hbar ^2}(p')^2} e^{-\mathrm {i}\frac{x' p'}{\hbar }}\nonumber \\&\quad \bigg \langle {x-\frac{1}{2}x'\bigg |{\hat{\rho }} \bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' p}{\hbar }}}\nonumber \\&= \frac{1}{2\pi \hbar }\sqrt{\frac{1}{\pi \alpha \hbar ^2}}\int {\hbox {d}x'}\int {\hbox {d}p'}e^{-\frac{1}{\alpha \hbar ^2}(p')^2}\nonumber \\&\quad \bigg \langle {x-\frac{1}{2}x'\bigg |{\hat{\rho }} \bigg |x+\frac{1}{2}x'\bigg \rangle e^{\mathrm {i}\frac{x' (p-p')}{\hbar }}}\nonumber \\&= \frac{1}{\sqrt{\pi \alpha \hbar ^2}}\int {\hbox {d}p'}e^{-\frac{1}{\alpha \hbar ^2}(p')^2}W(x,p-p'). \end{aligned}$$

(12)

Consequently, the spontaneous multiplication by a Gaussian function of the position mediated by the operator ${\hat{T}}({\hat{\rho }})$ in the Hilbert space formalism corresponds to a convolution with a Gaussian function of the momentum in the Wigner function formalism.

We can now combine Eqs. (8), (9), and  (12) to the GRW master equation for Wigner functions, given by

$$\begin{aligned} \partial _t W(x,p)= & {}\left \{H(x,p),W(x,p)\right\}_\star - \lambda \bigg (W(x,p)\nonumber \\&\quad - \frac{1}{\sqrt{\pi \alpha \hbar ^2}}\int {\hbox {d} p'}e^{-\frac{1}{\alpha \hbar ^2}(p')^2}W(x,p-p')\bigg ). \end{aligned}$$

(13)

To simplify Eq. (13), we make two approximations. First, we use a Kramers–Moyal expansion [96], which at second order leads to a Fokker–Planck equation [97]. For this purpose, we Taylor expand $W(x,p-p')$ around $p'=0$ and find

$$\begin{aligned} \partial _t W(x,p)= & {} \left\{H(x,p),W(x,p)\right\}_\star \nonumber \\&\quad +\sum _{n=1}^{\infty } \lambda \frac{(\alpha \hbar ^2)^n(2n-1)!!}{2^n (2n)!}\partial _p^{2n}W(x,p) \end{aligned}$$

(14)

where ! is the factorial and !! is the double factorial. When truncating at order $n=1$, we obtain from Eq. (14) the Fokker–Planck equation

$$\begin{aligned} \partial _t W(x,p) = \{H(x,p),W(x,p)\}_\star + D_p \partial _p^2 W(x,p) \end{aligned}$$

(15)

with the momentum diffusion coefficient $D_p = \lambda \alpha \hbar ^2/4$. Equations of the same form (but with a different physical content) have been obtained in the context of decoherence [98]. Here, one typically also has a friction term (see Sect. 9), Eq. (15) is then obtained when taking the limit of vanishing friction at fixed diffusion constant [20]. Moreover, equations similar to Eq. (15) have been obtained for the collapse-induced [47] and decoherence-induced [99] dynamics of quantum rotors.

Since Eq. (15) is of order $\hbar ^2$, the first term in Eq. (13) (the Moyal bracket) should also be expanded up to this order. This gives

$$\begin{aligned}\{H(x,p),W(x,p)\}_\star& \approx \nonumber - \frac{p}{m}\partial _x W(x,p) + (\partial _x U(x)) (\partial _p W(x,p))\nonumber \\&\quad - \frac{\hbar ^2}{24}(\partial _x^3 U(x)) (\partial _p^3 W(x,p)). \end{aligned}$$

(16)

We have assumed here that H has the form $H(x,p) = p^2/(2m) + U(x)$ with a potential U. If U is quadratic, linear, or constant, e.g., if we have a free particle or a harmonic oscillator, the time evolution obtained from Eq. (16) is thus classical. Making this assumption, the Fokker-Planck equation (15) reduces to

$$\begin{aligned}\partial _t W(x,p)& = - \frac{p}{m}\partial _x W(x,p) \nonumber \\&\quad + (\partial _x U(x)) (\partial _p W(x,p)) + D_p \partial _p^2 W(x,p). \end{aligned}$$

(17)

The same result is obtained if, e.g., due to a weak potential, the last term in Eq. (16) is small compared to the diffusion term in Eq. (17).

Now, we consider N particles in three dimensions. A calculation analogous to Eq. (12) gives

$$\begin{aligned} \sum _{i=1}^{N}T_i(W(\{\vec {r}_j,\vec {p}_j\}))= & {} \sum _{i=1}^{N} \frac{1}{(\pi \alpha \hbar ^2)^{\frac{3}{2}}}\nonumber \\&\quad \int {\hbox {d}^3 p'}e^{-\frac{1}{\alpha \hbar ^2}(\vec {p}')^2}\nonumber \\&\quad W(\{\vec {r}_j\},\vec {p}_1,\dotsc ,\vec {p}_{i-1},\vec {p}_i -\vec {p}',\vec {p}_{i+1},\dotsc ). \end{aligned}$$

(18)

Hence, the extension of Eq. (13) to N particles in three dimensions is given by

$$\begin{aligned}x\partial _t W(\{\vec {r}_j,\vec {p}_j\}) &=\left \{H(\{\vec {r}_j,\vec {p}_j\right\}), W(\{\vec {r}_j,\vec {p}_j\})\}_\star \nonumber \\&\quad -\lambda \bigg (W(\{\vec {r}_j,\vec {p}_j\}) - \sum _{i=1}^{N}T_i(W(\{\vec {r}_j,\vec {p}_j\}))\bigg ) \end{aligned}$$

(19)

with $\sum _i T_i$ given by Eq. (18). The Fokker–Planck equation (17) generalizes to

$$\begin{aligned}&\partial _t W(\{\vec {r}_j,\vec {p}_j\}) = \sum _{i=1}^{N}\nonumber \\&\quad \bigg (- \frac{\vec {p}_i}{m}\cdot \vec {\nabla }_{\vec {r}_i}W(\{\vec {r}_j,\vec {p}_j\}) + (\vec {\nabla }_{\vec {r}_i} U(\{\vec {r}_j\}))\cdot \nonumber \\&\quad (\vec {\nabla }_{\vec {p}_i} W(\{\vec {r}_j,\vec {p}_j\})) + D_p \vec {\nabla }_{\vec {p}_i}^2 W(\{\vec {r}_j,\vec {p}_j\})\bigg ) \end{aligned}$$

(20)

with $D_p = 3\lambda \alpha \hbar ^2/4$. This derivation shows that the Wigner approach is very useful in the study of spontaneous collapse models: Using Hilbert space theory, the derivation of a Fokker-Planck equation for GRW theory (as done in Ref. [38]) requires several pages of calculation and a significant number of auxiliary assumptions. In the Wigner formalism, however, it is obtained almost immediately and in a very natural way. On the other hand, it is notable that Eq. (13) has the form of a convolution in momentum space, whereas in Eq. (9) the density operator is simply multiplied by a Gaussian in position space. For this reason, if we do not operate with a simplified (Fokker-Planck) model, it is easier to work with the Fourier-transformed Wigner function

$$\begin{aligned} {\tilde{W}}(x,x') = \int {\hbox {d}p}W(x,p)e^{-\mathrm {i}\frac{x'p}{\hbar }} \end{aligned}$$

(21)

that depends on two position coordinates rather than on a position and a momentum coordinate. A physical reason for why it is easier to work with the Fourier-transformed Wigner function (21) is that Wigner functions typically treat position and momentum on an equal footing, whereas in GRW theory collapses take place in position space. For the Fourier-transformed Wigner function (21), Eq. (12) simplifies to

$$\begin{aligned} T({\tilde{W}}(x,x'))= e^{-\frac{\alpha }{4}(x')^2}{\tilde{W}}(x,x'). \end{aligned}$$

(22)

The result (22) also directly follows from Eq. (11) together with the fact that ${\tilde{W}}$ is related to the density operator in position representation, $\rho (x,x')=\langle x|{\hat{\rho }}|x'\rangle$, by ${\tilde{W}}(x,x')=\rho (x - x'/2, x+ x'/2)$ [as is obvious from Eq. (2)]. The relation of W and $\rho$ via the Fourier transformation [100] also explains why we need a convolution when describing the GRW dynamics in terms of Wigner functions (a Fourier transformation converts a multiplication to a convolution).

5 CSL model

A drawback of the GRW theory is that it does not preserve the symmetry properties of a system of identical particles [90]. In particular, it does not respect the fact that the wavefunction of a system of identical particles has to be symmetric under particle exchange for bosons and antisymmetric for fermions. This problem was solved by the continuous spontaneous localization (CSL) model introduced in Refs. [39, 63]. We here follow Ref. [90]. (A discussion of this model can also be found in Refs. [101, 102].) First, we define the locally averaged density operator

$$\begin{aligned} {\hat{N}}(\vec {r})=\sum _{s}^{}\int {\hbox {d}^3r'}g(\vec {r}'-\vec {r}){\hat{a}}^\dagger (\vec {r}',s){\hat{a}} (\vec {r}',s), \end{aligned}$$

(23)

where ${\hat{a}}^\dagger (\vec {r}',s)$ and ${\hat{a}}^\dagger (\vec {r}',s)$ are creation and annihilation operators for a particle with spin component s at position $\vec {r}'$ and g is a symmetric positive function typically chosen to be

$$\begin{aligned} g(\vec {r})=\Big (\frac{\alpha }{2\pi }\Big )^\frac{3}{2}\exp \!\Big (-\frac{\alpha }{2} \vec {r}^2\Big ). \end{aligned}$$

(24)

The density operator is assumed to follow the dynamic equation

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}{\hat{\rho }}&= - \frac{\mathrm {i}}{\hbar }[{\hat{H}},{\hat{\rho }}] \nonumber \\&\quad \,+\xi \int {\hbox {d}^3r}\Big ({\hat{N}}(\vec {r}){\hat{\rho }}{\hat{N}}(\vec {r}) -\frac{1}{2}[ {\hat{N}}^2(\vec {r}),{\hat{\rho }}]_+\Big ) \end{aligned}$$

(25)

with a constant $\xi$ that is related to the localization rate, see Eq. (29) below, and the anticommutator $[\cdot ,\cdot ]_+$. From here on, we ignore the spin degrees of freedom since they are not relevant for the collapse dynamics. See Ref. [90] for a discussion of how Eq. (25) follows from the theory of stochastic processes in Hilbert space.

We now derive the Wigner representation of Eq. (25). In position representation, the last two terms on the right-hand side of Eq. (25) read [90]

$$\begin{aligned}&\frac{\xi }{2}\bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}} \sum _{i,j=1}^{N}\langle \vec {r}_1,\dotsc ,\vec {r}_N|{\hat{\rho }}| \vec {r}_1',\dotsc ,\vec {r}_N'\rangle \nonumber \\&\quad \Big (2e^{-\frac{\alpha }{4}(\vec {r}_i-\vec {r}_j')^2} -e^{-\frac{\alpha }{4}(\vec {r}_i-\vec {r}_j)^2}-e^{-\frac{\alpha }{4}(\vec {r}_i'-\vec {r}_j')^2}\Big ). \end{aligned}$$

(26)

Thus, the Weyl symbol for the operator given by the last two terms on the right-hand side of Eq. (25) is

$$\begin{aligned}&\frac{1}{(2\pi \hbar )^{3N}}\int {\hbox {d}^3r_1'}\cdots \int {\hbox {d}^3r_N'} \frac{\xi }{2}\bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}}\nonumber \\&\quad \sum _{i,j=1}^{N} \bigg \langle \vec {r}_1-\frac{\vec {r}_1'}{2},\dotsc ,\vec {r}_N-\frac{\vec {r}_N'}{2} \bigg |{\hat{\rho }}\bigg |\vec {r}_1 +\frac{\vec {r}_1}{2},\dotsc ,\vec {r}_N+\frac{\vec {r}_N'}{2}\bigg \rangle \nonumber \\&\,e^{\frac{\mathrm {i}}{\hbar }\sum _{k=1}^{N}\vec {r}_k'\cdot \vec {p}_k}\nonumber \\&\quad \Big (2e^{-\frac{\alpha }{4}(\vec {r}_i-\frac{1}{2}\vec {r}_i'-\vec {r}_j-\frac{1}{2}\vec {r}_j')^2} -e^{-\frac{\alpha }{4}(\vec {r}_i-\frac{1}{2}\vec {r}_i'-\vec {r}_j+\frac{1}{2}\vec {r}_j')^2}\nonumber \\&\quad -e^{-\frac{\alpha }{4}(\vec {r}_i+\frac{1}{2}\vec {r}_i'-\vec {r}_j-\frac{1}{2}\vec {r}_j')^2}\Big )\nonumber \\&=\frac{1}{(2\pi \hbar )^{3N}}\bigg (\frac{1}{\pi \alpha \hbar ^2}\bigg )^\frac{3}{2}\int {\hbox {d}^3r_1'}\cdots \int {\hbox {d}^3r_N'}\int {\hbox {d}^3 p'}\frac{\xi }{2}\bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}}\nonumber \\&\quad \sum _{i,j=1}^{N}\bigg \langle \vec {r}_1-\frac{\vec {r}_1'}{2},\dotsc ,\vec {r}_N-\frac{\vec {r}_N'}{2} \bigg |{\hat{\rho }}\bigg |\vec {r}_1 +\frac{\vec {r}_1}{2},\dotsc ,\vec {r}_N+\frac{\vec {r}_N'}{2}\bigg \rangle \nonumber \\&\quad \; e^{\frac{\mathrm {i}}{\hbar }\sum _{k=1}^{N}\vec {r}_k'\cdot \vec {p}_k}e^{-\frac{(\vec {p}')^2}{\alpha \hbar ^2}}\nonumber \\&\quad \Big (2e^{\frac{\mathrm {i}}{\hbar }(\vec {r}_i-\frac{1}{2}\vec {r}_i'-\vec {r}_j-\frac{1}{2}\vec {r}_j')\cdot \vec {p}'} -e^{\frac{\mathrm {i}}{\hbar }(\vec {r}_i-\frac{1}{2}\vec {r}_i'-\vec {r}_j+\frac{1}{2}\vec {r}_j')\cdot \vec {p}'}\nonumber \\&\quad -e^{\frac{\mathrm {i}}{\hbar }(\vec {r}_i+\frac{1}{2}\vec {r}_i'-\vec {r}_j-\frac{1}{2}\vec {r}_j')\cdot \vec {p}'}\Big )\nonumber \\&=\bigg (\frac{1}{\pi \alpha \hbar ^2}\bigg )^\frac{3}{2}\int {\hbox {d}^3 p'} \frac{\xi }{2}\bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}} e^{-\frac{(\vec {p}')^2}{\alpha \hbar ^2}}\nonumber \\&\quad \sum _{i,j=1}^{N}e^{\frac{\mathrm {i}}{\hbar }(\vec {r}_i-\vec {r}_j)\cdot \vec {p}'}\nonumber \\&\quad \bigg (2W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i -\frac{\vec {p}'}{2},\vec {p}_j-\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg ) \nonumber \\&\quad \,-W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i-\frac{\vec {p}'}{2},\vec {p}_j+\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg )\nonumber \\&\quad -W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i+\frac{\vec {p}'}{2},\vec {p}_j-\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg )\!\bigg ). \end{aligned}$$

(27)

Therefore, the CSL master equation for Wigner functions reads

$$\begin{aligned} \partial _t W(\{\vec {r}_i,\vec {p}_i)\}&=\left \{H(\{\vec {r}_i,\vec {p}_i\}),W(\{\vec {r}_i,\vec {p}_i\}\right\}_\star \nonumber \\&\quad \,+\bigg (\frac{1}{\pi \alpha \hbar ^2}\bigg )^\frac{3}{2}\int {\hbox {d}^3 p'}\frac{\xi }{2}\bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}}e^{-\frac{(\vec {p}')^2}{\alpha \hbar ^2}}\nonumber \\&\quad \sum _{i,j=1}^{N}e^{\frac{\mathrm {i}}{\hbar }(\vec {r}_i-\vec {r}_j)\cdot \vec {p}'}\nonumber \\&\quad \bigg (2W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i-\frac{\vec {p}'}{2},\vec {p}_j-\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg ) \nonumber \\&\quad \,-W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i-\frac{\vec {p}'}{2},\vec {p}_j+\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg )\nonumber \\&\quad -W\bigg (\{\vec {r}_k\},\vec {p}_1,\dotsc ,\vec {p}_i+\frac{\vec {p}'}{2},\vec {p}_j-\frac{\vec {p}'}{2},\dotsc ,\vec {p}_N\bigg )\!\bigg ). \end{aligned}$$

(28)

Equations (28) and (19) coincide for $N=1$ if we set [90]

$$\begin{aligned} \lambda = \xi \bigg (\frac{\alpha }{4\pi }\bigg )^{\frac{3}{2}}. \end{aligned}$$

(29)

6 Quantum gravity

We now come back to the measurement problem. An interesting solution in the spirit of dynamical collapse models was proposed by Penrose [64, 65] and Diósi [66, 67], which is based on the problem of quantum gravity: As is well known, there is at present no generally accepted and experimentally verified unification of quantum mechanics and general relativity. While most approaches to this issue are based on quantizing gravity, an alternative approach consists in “gravitizing quantum mechanics” [64], i.e., in modifying quantum mechanics to incorporate effects of gravity. A spatial quantum superposition generates a superposition of spacetimes, which is supposed to lead to a gravity-induced collapse of the wavefunction [45]. Since this effect will depend on the mass of the system, it will be important for macroscopic, but negligible for microscopic systems; explaining why superpositions are not observed on macroscopic scales. Recent experimental tests [45], however, found no evidence of a gravity-induced collapse of the wave function, putting significant constraints on the possible validity of the Diósi-Penrose model.

A simple quantitative model for this effect is given by [45, 66, 67]

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}{\hat{\rho }} = - \frac{\mathrm {i}}{\hbar }[{\hat{H}},{\hat{\rho }}] - \frac{4\pi G}{\hbar }\int {\hbox {d}^3r}\int {\hbox {d}^3r'} \frac{[{\hat{M}}(\vec {r}'),[{\hat{M}}(\vec {r}),{\hat{\rho }}]]}{\Vert \vec {r}-\vec {r}'\Vert } \end{aligned}$$

(30)

with the gravitational constant G, the total mass density operator ${\hat{M}}$, and the Euclidean norm $\Vert \cdot \Vert$. This operator is given by [45]

$$\begin{aligned} {\hat{M}}(\vec {r})=\sum _{i=1}^{N}{\hat{\mu }}_i(\vec {r}) \end{aligned}$$

(31)

with the mass-density operator ${\hat{\mu }}_i$ for the i-th particle. The Diósi-Penrose model constitutes a third type of collapse models in addition to GRW theory and CSL model.

A difficulty arises from the fact that the standard definition

$$\begin{aligned} {\hat{\mu }}_i(\vec {r})=m_i\delta (\vec {r}-\hat{\vec {r}}_i) \end{aligned}$$

(32)

with the mass $m_i$ of and the position operator $\hat{\vec {r}}_i$ for the i-th particle leads to divergences [45]. To solve this problem, one has to smear out the mass density. This can be done in various ways [93]. Proposals from the literature include [67]

$$\begin{aligned} {\hat{\mu }}_i(\vec {r})=\frac{3m_i}{4\pi R_0^2}\Theta (R_0 - \Vert \vec {r}-\hat{\vec {r}}_{i}\Vert ), \end{aligned}$$

(33)

with a length scale $R_0$ and the Heaviside function $\Theta$, as well as [103]

$$\begin{aligned} {\hat{\mu }}_i(\vec {r})=\frac{m_i}{(2\pi R_0^2)^{\frac{3}{2}}}\exp \!\bigg (-\frac{(\vec {r}-\hat{\vec {r}}_i)^2}{2R_0^2}\bigg ). \end{aligned}$$

(34)

We are only interested in the non-Hamiltonian term on the right-hand side of Eq. (30). Moreover, we restrict ourselves to a single particle with mass density operator ${\hat{\mu }}$. As shown in the supplementary material of Ref. [45], the non-Hamiltonian term in Eq. (30) can be written as

$$\begin{aligned} \int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')\Big (e^{\frac{\mathrm {i}}{\hbar }\vec {p}'\cdot \hat{\vec {r}}}{\hat{\rho }}\, e^{-\frac{\mathrm {i}}{\hbar }\vec {p}'\cdot \hat{\vec {r}}} - {\hat{\rho }}\Big ) \end{aligned}$$

(35)

with the function

$$\begin{aligned} {\tilde{\Gamma }}(\vec {p}) =\frac{4 G}{\pi \hbar ^2}\frac{|{\tilde{\mu }}(\vec {p})|^2}{\Vert \vec {p}\Vert ^2} \end{aligned}$$

(36)

depending on the Fourier-transformed mass density ${\tilde{\mu }}$. The matrix element of Eq. (35) is given by [45]

$$\begin{aligned} \int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')\Big (e^{\frac{\mathrm {i}}{\hbar }\vec {p}'\cdot (\vec {r}-\vec {r}')}-1\Big ) \langle \vec {r}|{\hat{\rho }}|\vec {r}'\rangle . \end{aligned}$$

(37)

Hence, the Weyl symbol of the expression (35) is given by

$$\begin{aligned}&\frac{1}{(2\pi \hbar )^3}\int {\hbox {d}^3r'}\int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')\Big (e^{-\frac{\mathrm {i}}{\hbar }\vec {p}'\cdot \vec {r}'}-1\Big )\nonumber \\&\bigg \langle \vec {r}-\frac{\vec {r}'}{2} \bigg |{\hat{\rho }}\bigg |\vec {r} + \frac{\vec {r}'}{2}\bigg \rangle e^{\frac{\mathrm {i}}{\hbar }\vec {p}\cdot \vec {r}'}\nonumber \\&=\frac{1}{(2\pi \hbar )^3}\int {\hbox {d}^3r'}\int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')\bigg \langle \vec {r} -\frac{\vec {r}'}{2}\bigg |{\hat{\rho }}\bigg |\vec {r} + \frac{\vec {r}'}{2}\bigg \rangle e^{\frac{\mathrm {i}}{\hbar }(\vec {p}-\vec {p}')\cdot \vec {r}'} \nonumber \\&\quad \, - W(\vec {r},\vec {p})\nonumber \\&=\int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')(W(\vec {r},\vec {p}-\vec {p}') - W(\vec {r},\vec {p})). \end{aligned}$$

(38)

Using Eq. (38), the Diósi-Penrose master equation for Wigner functions is given by

$$\begin{aligned} \partial _t W(\vec {r},\vec {p})&= \left\{H(\vec {r},\vec {p}), W(\vec {r},\vec {p})\right\}_\star \nonumber \\&\quad \, + \int {\hbox {d}^3p'}{\tilde{\Gamma }}(\vec {p}')(W(\vec {r},\vec {p}-\vec {p}') - W(\vec {r},\vec {p})). \end{aligned}$$

(39)

When comparing Eq. (39) with Eq. (13), it is easily seen that the general structure is the same, in particular if we assume that the mass-density distribution is Gaussian. However, the expression that the Wigner function is convoluted with differs due to the presence of the factor $1/(\Vert \vec {p}\Vert ')^2$ in Eq. (36). A similar observation was made in Ref. [93]. Equation (39) is an important result, since it has recently been suggested that Wigner functions could be helpful for experimental tests of quantum gravity [104].

7 Dissipative GRW model

Finally, we turn to a more recent extension of the GRW theory, the dissipative GRW model proposed by Smirne et al. [68]. As discussed in Sect. 4, the standard spontaneous collapse models have the problem that the energy increases continuously. This increase is not measurable on typical observational timescales, but leads to a divergence of the temperature for $t\rightarrow \infty$. To avoid this problem, the dissipative GRW model includes dissipation to ensure that the system thermalizes to a finite temperature. A similar extension has been developed for the CSL model by Smirne and Bassi [69] and for the Diósi-Penrose model by Bahrami et al. [93].

Since the derivation of the Wigner representation is more involved for this theory, we restrict the discussion to a single particle in one dimension to improve readability (the extension is straightforward). The dissipative GRW model is given by the master equation [68]

$$\begin{aligned} \partial _t {\hat{\rho }} = \frac{\mathrm {i}}{\hbar }[{\hat{H}},{\hat{\rho }}] - \lambda ({\hat{\rho }} - {\hat{T}}_\mathrm {d}({\hat{\rho }})), \end{aligned}$$

(40)

which looks exactly like Eq. (9). However, the standard collapse operator ${\hat{T}}$ is replaced by the dissipative collapse operator

$$\begin{aligned}&{\hat{T}}_\mathrm {d}({\hat{\rho }}) = \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\int {\hbox {d}p'}e^{\frac{\mathrm {i}}{\hbar }p'{\hat{x}}}e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}\nonumber \\&\quad {\hat{\rho }}\, e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}e^{-\frac{\mathrm {i}}{\hbar }p'{\hat{x}}} \end{aligned}$$

(41)

with a localization length scale $r_{\mathrm {c}} = 1/\sqrt{\alpha }$ and a temperature parameter k. The operator (41) reduces to the standard collapse operator (10) (and the dissipative GRW model [40) to the ordinary GRW model (9)] for $k \rightarrow 0$, which corresponds to the high-temperature limit.

A transformation of Eqs. (40) and (41) to the Wigner formalism is carried out in “Appendix A”. The result, the dissipative GRW master equation for Wigner functions, reads

$$\begin{aligned}\partial _t W(x,p) &= \{H(x,p),W(x,p)\}_\star \nonumber \\&\quad - \lambda (W(x,p) - T_\mathrm {d}(W(x,p))) \end{aligned}$$

(42)

with

$$\begin{aligned}&T_\mathrm {d}(W(x,p))= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{\sqrt{2\pi k^2 r_{\mathrm {c}}^2}}\nonumber \\&\quad \int {\hbox {d}p'}\int {\hbox {d}x'} e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8k^2p^2 + 8 k p p' - 8 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2)}\nonumber \\&\quad e^{-\frac{(x')^2}{2k^2 r_{\mathrm {c}}^2}}W(x-x',p-p'). \end{aligned}$$

(43)

Performing a Kramers–Moyal expansion in analogy to the treatment in Sect. 4 gives (as shown in “Appendix A”) the Fokker–Planck equation

$$\begin{aligned}\partial _t W(x,p)& = \{H(x,p),W(x,p)\}_\star + \gamma \partial _p (p W(x,p))\nonumber \\&\quad + \gamma m k_{\mathrm {B}} T_\mathrm {n} \partial _p^2 W(x,p) \end{aligned}$$

(44)

with damping parameter $\gamma = 2\lambda k$, Boltzmann constant $k_{\mathrm {B}}$, and noise temperature $T_\mathrm {n} = \hbar ^2 /(8 m k k_{\mathrm {B}} r_{\mathrm {c}}^2)$. This is precisely the temperature a dissipative GRW system thermalizes to [105], which confirms the consistency of our approach. If we ignore the fact that Eq. (44) contains a Moyal rather than a Poisson bracket, it is formally identical to the Kramers equation [106]

$$\begin{aligned} \partial _t W(x,p)= & {} \Big (- \frac{p}{m}\partial _x + (\partial _x U(x)) \partial _p + \partial _p\gamma p + \gamma m k_{\mathrm {B}} T_\mathrm {b} \partial _p^2 \Big )\nonumber \\&\quad W(x,p), \end{aligned}$$

(45)

which describes the dynamics of an underdamped Brownian particle in contact with a heat bath that has temperature $T_\mathrm {b}$. In Eq. (44), however, the damping and diffusion terms arise not from the interaction with a solvent, but from the fundamental equations of motion of (GRW-type) quantum mechanics, implying that they are also present in a closed system.

8 Irreversibility and approach to equilibrium

We now move to a different topic, namely the foundations of statistical mechanics. Among the central problems in this field is the question how, given the fact that the microscopic laws of physics are time-reversal invariant, macroscopic thermodynamic irreversibility can be explained [52, 53]. A widely discussed idea in this context is David Albert’s [54,55,56] suggestion that the spontaneous collapses postulated by GRW theory might play a role in the explanation of irreversibility. In the remainder of this article, we make use of the results from Sects. 4 and 7 to discuss this proposal and to test (and refute) it via numerical simulations. The significance of this test lies in two aspects. First, it constitutes an application of the results from the previous sections and thereby demonstrates the usefulness of the Wigner approach for GRW theory. Second, and more importantly, it constitutes the first quantitative examination of a well-known proposal from philosophy of physics. The null result presented here has to be taken into account when discussing GRW-based approaches to irreversibility.

We now introduce Albert’s proposal following Refs. [54,55,56]. A typical explanation for the approach to thermodynamic equilibrium runs as follows: Macroscopic states of a system are associated with a large number of possible microscopic states. The equilibrium state is the macrostate with the largest entropy, i.e., the state with the largest number of possible microscopic realizations. Consequently, most of the microstates associated with an initial macrostate will evolve towards an equilibrium state, and since we do not know the exact microstate of the system, we have to conclude that it is considerably more likely that the system will evolve towards equilibrium rather than away from it.

This explanation, Albert argues, is unsatisfactory because it involves an appeal to the limitations of observers (ignorance about the precise initial microstate) in explaining what appears to be an observer-independent phenomenon, namely the approach to equilibrium. GRW theory provides a possible solution by introducing objective probabilities associated with the spontaneous collapse of the wavefunction. If the system is initially in an “abnormal” microstate that would lead to an evolution away from equilibrium—which is unlikely, but not impossible—then a GRW collapse will in most cases lead to a “normal” state from which the system equilibrates. According to Albert, this happens since every microscopic neighborhood of the abnormal microstate in phase space contains significantly more normal than abnormal microstates, such that thermodynamic abnormality is unstable under perturbations. Some discussions of Albert’s proposal can be found in Refs. [57,58,59]. See Ref. [62] for a review.

Albert’s approach belongs to the so-called Boltzmann approach to statistical mechanics [58]. In this framework, one is interested in a single system that has a certain macrostate (characterized by its macroscopic properties). This macrostate depends on the microscopic state of the system. Various microstates can correspond to the same macrostate. The equilibrium state is then the macrostate that the largest number of microstates correspond to. An alternative framework, which is more widely used in practical applications, is the Gibbs approach, where one studies ensembles (hypothetical sets of infinitely many copies of a system with different initial conditions). “Approach to equilibrium” then means that the probability distribution over an ensemble of systems converges towards the equilibrium (in the simplest case: Maxwell-Boltzmann) distribution. A more detailed discussion and comparison of both approaches can be found in Ref. [107].

Albert’s proposal is conceptually appealing for at least three reasons [58]: First, it allows to solve two problems—the quantum measurement problem and the problem of thermodynamic irreversibility—at once. Second, it unifies quantum-mechanical probabilities with the probabilities used in statistical mechanics. Third, it allows to explain thermodynamic irreversibility without requiring assumptions about a probability distribution over initial conditions. However, it is only based on a qualitative argument, a quantitative test is lacking at present. Such a test would be required to see how far GRW collapses can actually get us on the road to equilibrium [62].

A simple argument for Albert’s proposal can be developed based on Eq. (20). As is well known, the Boltzmann equation [108] provides a useful model for the approach to thermodynamic equilibrium. The one-particle distribution function $f(\vec {p})$, giving the probability that a particle has momentum $\vec {p}$, approaches the Maxwell-Boltzmann (equilibrium) distribution due to the presence of collisions. This is a consequence of the fact that, in the derivation of the Boltzmann equation, one assumes that velocities are uncorrelated, i.e., that the two-particle distribution satisfies $f_2(\vec {p}_1,\vec {p}_2) = f(\vec {p}_1)f(\vec {p}_2)$ (molecular chaos), before a collision [109]. If one instead assumes that velocities are uncorrelated after a collision, we would have found that the entropy always decreases.

Thus, if Eq. (20) should explain the approach to thermodynamic equilibrium, it should enforce $f_2(\vec {p}_1,\vec {p}_2) = f(\vec {p}_1)f(\vec {p}_2)$. There is a good reason to believe that this is the case: The distributions f and $f_2$ can be obtained by integrating the Wigner function³ over all positional and over all but one (for f) or two (for $f_2$) momentum degrees of freedom [21]. Ignoring the reversible Hamiltonian part of the dynamics, we find that the difference between a factorized and a nonfactorized distribution evolves as

$$\begin{aligned} \begin{aligned}&\partial _t (f_2(\vec {p}_1,\vec {p}_2) - f(\vec {p}_1)f(\vec {p}_2))\\&=D_p(\vec {\nabla }_{\vec {p}_1}^2 + \vec {\nabla }_{\vec {p}_2}^2)(f_2(\vec {p}_1,\vec {p}_2) - f(\vec {p}_1)f(\vec {p}_2)). \end{aligned} \end{aligned}$$

(46)

Equation (46) is a diffusion equation for the function ${\mathfrak {f}}(\vec {p}_1,\vec {p}_2)=f_2(\vec {p}_1,\vec {p}_2) - f(\vec {p}_1)f(\vec {p}_2)$. Due to the diffusive dynamics, ${\mathfrak {f}}$ will be homogeneous at long times, implying $f_2(\vec {p}_1,\vec {p}_2) = f(\vec {p}_1)f(\vec {p}_2) +c$ with a constant c. Since f and $f_2$ are normalized, we have $c=0$. Consequently, the GRW dynamics indeed tends to destroy correlations in momentum space and thereby leads to molecular chaos.

There is also a different aspect that has to be taken into account: When comparing Eq. (20) to the Kramers equation (45) we can note that Eq. (20) has no friction term. However, the friction term in Eq. (45) is known to be essential for the approach to an equilibrium state (corresponding, in the classical case, to the Maxwell-Boltzmann distribution), and is intimately connected to the fluctuations described by the diffusion term in Eq. (45) via a fluctuation-dissipation theorem. The fluctuations in Eq. (20), in contrast, are not connected to any form of dissipation. Thus, if Eq. (20) induces equilibration, it will be of a different form than equilibration processes resulting from friction.

This brings us back to the distinction between Boltzmannian and Gibbsian approaches to equilibrium introduced above. If W can be interpreted as a probability distribution (i.e., in the classical case), the Kramers equation (45) describes an approach towards the Maxwell-Boltzmann distribution, which is the equilibrium distribution of Gibbsian statistical mechanics. For this, the friction term plays an essential role. In Eq. (20), in contrast, a friction term is not present, such that we cannot expect that W converges to a Maxwell-Boltzmann distribution. This, however, does not exclude that it leads (with high probability) to a convergence towards the macrostate with the largest phase-space volume, i.e., to Boltzmannian equilibrium.

Up to now, we have only considered the original GRW model by Ghirardi et al. [38]. This is the theory that Albert’s original suggestion is based on.⁴ The drawbacks of the standard GRW model in explaining equilibration might, however, be avoided by considering the dissipative GRW model introduced in Sect. 7. As shown there, the Wigner function approximately follows the dynamic Eq. (44) in this theory. The fact that Eq. (44) does contain a damping term related to the diffusion term in the standard way (actually, Eq. (44) reduces to Eq. (45) for $\hbar \rightarrow 0$ if we set $T_\mathrm {n} = T_\mathrm {b}$) seems to suggest a modification of Albert’s proposal in which the approach to thermodynamic equilibrium is described by the dissipative GRW model.

This would have interesting consequences: The Kramers equation (45), of which Eq. (44) is an extension, describes a system coupled to a heat bath at fixed temperature (rather than an isolated system with fixed energy) that should be described in the canonical ensemble. Since spontaneous GRW-type collapses occur also in isolated systems, this would imply that the equilibrium distribution of an isolated system is given by the canonical rather than by the microcanonical distribution [69, 93]. This is a very significant difference to “standard” statistical mechanics, where isolated systems are described by a microcanonical distribution, which corresponds to the energy-conserving case. While the different ensembles agree in the thermodynamic limit, they can differ significantly for systems with a small number of particles [111]. A different result would be obtained in an energy-conserving spontaneous collapse model that might be derived in the future. (Note, however, that Albert’s original discussion [54, 55] is set up in a Boltzmannian framework that is not based on ensembles at all.) In addition, the reduction of the spontaneous heating by dissipation has implications for tests of spontaneous collapse models that rely on this effect, for example those based on measurements of neutron star heating [112].

The idea that the dissipative GRW model (44) is more appropriate as an explanation for thermodynamic irreversibility than the original GRW model is, however, solely based on the fact that it is formally analogous to Eq. (45), which is known to lead to equilibration. Physically, the dissipative GRW model essentially assumes that all particles in the universe are coupled to a universal heat bath with a certain temperature $T_\mathrm {n}$. (The nondissipative GRW model is the special case $T_\mathrm {n} \rightarrow \infty$.) Consequently, in the dissipative GRW model, systems will (in the long run) approach a canonical equilibrium state with temperature $T_\mathrm {n}$. This is irreversibility, but not of the type known from classical thermodynamics. There, a typical manifestation of equilibration is that if two systems with initially different temperatures $T_1$ and $T_2$ are brought in thermal contact, heat will flow from the hotter to the colder system until they have reached equal temperatures [113]. This final temperature will depend on the initial temperatures and other properties of the systems. It is certainly not a prediction of classical thermodynamics that all isolated systems converge to the same equilibrium state characterized by a universal temperature $T_\mathrm {n}$.

This, of course, does not necessarily imply that the dissipative GRW model is incorrect. It merely implies that it makes a different empirical prediction than standard thermodynamics, namely that there is a universal temperature that all systems, irrespective of their initial temperature and energy, will equilibrate to. This equilibration will take place on a very long timescale, such that this prediction does not contradict the observation of “normal” equilibration in the form of heat exchange on ordinary timescales. However, it was this equilibration that we were seeking an explanation for. In this regard, the dissipative GRW model will play the same role (if any) as the standard one, namely that of erasing correlations by providing dissipation.

These conclusions are not affected in a significant way by working (as Albert does) in a Boltzmannian rather than a Gibbsian framework. As discussed by te Vrugt [53], the distinction between Gibbsian and Boltzmannian statistical mechanics can be linked to the distinction between deterministic and stochastic forms of dynamical density functional theory (DDFT) [33]. DDFT is a theory for the nonequilibrium evolution of the one-body density in classical fluids. Deterministic DDFT describes ensembles, whereas stochastic DDFT a (possibly spatially averaged) single system [114]. We can thus get an idea of the Boltzmannian equilibration predicted by dissipative GRW models by considering the predictions of stochastic DDFT, namely an approach to the minimum state of the free energy functional (which has a different definition compared to deterministic DDFT [33]). Hence, according to dissipative GRW theory, the system will (in equilibrium) fluctuate around the minimum of a free energy functional⁵ with temperature $T_\mathrm {n}$.

When it comes to solving the fundamental problems of explaining thermodynamic irreversibility, typical arguments against coarse-graining—namely that the equilibrium state of equilibrium statistical mechanics is reproduced only approximately [117] and on certain timescales—also apply to GRW-based explanations. Nevertheless, Albert’s approach, if successful, would still explain why anti-thermodynamic behavior is never observed in macroscopic systems. For microscopic systems, it has recently been demonstrated experimentally that heat can spontaneously flow from a cold to a hot system for initially correlated spins [118]. This is not in contradiction to the GRW-based approach, since GRW collapses are (by construction) extremely rare for small quantum systems and thus allow for such a behavior.

9 Connection to decoherence

Spontaneous collapse models like the GRW theory have a very close connection to the theory of quantum decoherence [119], which is reviewed in Refs. [77, 79, 98, 120]. In fact, some authors refer to GRW-type collapses as “intrinsic decoherence” [121]. Moreover, there has been a significant amount of research on explaining the emergence of irreversibility based on the interaction of quantum systems with their environment [122], in particular decoherence [77, 120, 123]. Consequently, it is of interest to relate such approaches to Albert’s proposal discussed here. The purpose of this section (which is a slight digression) is thus to present some other approaches from the literature that Albert’s theory, on which we focus here, is related to. This allows to see our results in a more general context.

In contrast to the GRW theory, decoherence does by itself neither provide a solution to the measurement problem nor produce stochasticity. However, it is an important ingredient in various interpretations of quantum mechanics [81]. This has been exploited in an explanation of the approach to equilibrium proposed by Hemmo and Shenker [57,58,59]. Like Albert’s idea, it is based on explaining probabilities in statistical mechanics based on quantum-mechanical probabilities. However, Hemmo and Shenker take as a starting point not GRW theory, but no-collapse interpretations (such as modal interpretations or the Everett interpretation), in which a collapse of the wavefunction does not occur. Since the quantum dynamics of closed systems is isolated in the absence of collapses, stochasticity then requires external interventions. Thus, the proposal discussed in this section belongs to the tradition of interventionism [117, 124, 125], which explains the approach to equilibrium based on external perturbations.

We discuss the approach following Ref. [58]. In decoherence theory, one typically assumes that the initial state $|\Psi \rangle$ of a quantum system and its environment can be written in the form

$$\begin{aligned} |\Psi \rangle (0)=|\psi \rangle \otimes |E\rangle \end{aligned}$$

(47)

with the state of the system $|\psi \rangle$, the state of the environment $|E\rangle$, and the tensor product $\otimes$. Moreover, one assumes that the interaction Hamiltonian describing the interaction of the system with its environment commutes with a certain observable, the so-called pointer variable. A typical pointer variable is the position. Let $\{|\psi _i\rangle \}$ be a basis of the Hilbert space of the quantum system that consists of eigenstates of the pointer variable. The state of the system at time t can then be written as

$$\begin{aligned} |\Psi \rangle (t)=\sum _{i}\nu _i(t)|\psi _i\rangle \otimes |E_i\rangle \end{aligned}$$

(48)

with the relative states $|E_i\rangle$ and the amplitudes $\nu _i(t)$. Due to the coupling with the environment, the relative states will satisfy $\langle E_i|E_j\rangle \approx \delta _{ij}$ with the Kronecker delta $\delta _{ij}$ after an extremely short time. Then, the reduced density operator (obtained by tracing over environmental degrees of freedom) of the system takes the form

$$\begin{aligned} {\hat{\rho }}(t) = \sum _{i}|\nu _i(t)|^2 |\psi _i\rangle \langle \psi _i|. \end{aligned}$$

(49)

For a harmonic oscillator weakly coupled to an environment in equilibrium, the states $|\psi _i\rangle$ correspond to narrowly peaked Gaussians, such that the decohered system follows quasi-classical trajectories [126].

This result is, by itself, no solution of the quantum measurement problem since it leads to a superposition of different classical time evolutions. However, in no-collapse interpretations of quantum mechanics, such as modal interpretations [127], it is possible to interpret the reduced state (49) as a probability distribution over the different quasi-classical states $\{|\psi _i\rangle \}$. During the time evolution, there will then be stochastic transitions between such states. These stochastic transitions can then play the same role in the approach to equilibrium as the stochastic GRW collapses in Albert’s theory.

What is interesting about this approach is that the Wigner function of a decohering system coupled to a thermal environment can, in the simplest case, be shown to have the form [20, 58, 123, 128,129,130]

$$\begin{aligned} \begin{aligned} \partial _t W(x,p)&= \{H(x,p),W(x,p)\}_\star \\&\quad \, +(\partial _p\gamma p+ \gamma m k_{\mathrm {B}} T_\mathrm {b} \partial _p^2) W(x,p), \end{aligned} \end{aligned}$$

(50)

which, for $\hbar \rightarrow 0$, reduces to the Kramers equation (45). Thus, in this approach, the probability distribution also approaches a Maxwell-Boltzmann form in the classical case. Consequently, it is, like the dissipative GRW model, close in form to standard theories of equilibration. However, for these models, the temperature the system equilibrates to is that of the environment (rather than that of a universal heat bath).

It should be noted that there are important conceptual differences between the GRW-based and the decoherence-based explanation of equilibration—if we explain thermodynamic irreversibility by decoherence, then there can be no irreversibility in an isolated system. In contrast, such a behavior is to be expected from the GRW perspective, which would thus (in a sense) imply that the standard models of stochastic dynamics (such as the Fokker-Planck equation or the Langevin equation) are actually more fundamental than the ordinary Hamiltonian dynamics which they are usually thought to be an approximation to.

A further recent approach to the problem of explaining the emergence of statistical mechanics that is interesting in this context was proposed by Drossel [131,132,133,134]. She argues that the practice of statistical mechanics, including (but not limited to) the presence of irreversibility, is incompatible with the universal validity of deterministic unitary quantum mechanics, and has to be explained by the presence of fundamental stochasticity. Her view can be classified as an intermediate position between that of Albert and that of Hemmo and Shenker. Like Albert (and unlike Hemmo and Shenker), she sees the origin of irreversibility in a stochastic quantum dynamics that violates the standard Schrödinger time evolution. However, like Hemmo and Shenker (and unlike Albert), she attributes this stochasticity to interactions of a quantum system with the environment. Finally, Wallace [29] has argued, within an Everettian framework, that the probabilities of statistical mechanics can be understood as arising from quantum-mechanical probabilities based on the fact that the classical limit of quantum mechanics is approximately isomorphic to a theory of classical probability distributions (as is evident from the Wigner function picture).

10 Numerical experiments

Fig. 1

a Initial state of the first simulation (reference state). The system consists of two parts in thermal contact, separated by a sharp temperature gradient. b Final (relatively homogeneous) temperature distribution of the first simulation, obtained from evolving the system shown in (a) forward in time using Hamiltonian dynamics. c Final temperature distribution of the second simulation, obtained from evolving a system with “abnormal” initial conditions (state shown in (b) with reversed momenta) under Hamiltonian dynamics. The system moves away from thermal equilibrium. d Final temperature distribution of the third simulation, obtained from evolving a system with the same initial condition as in (c) using GRW dynamics. The stochastic perturbations make no difference for the anti-thermodynamic behavior. e Final temperature distribution of the fourth simulation, obtained from evolving a system with the same initial condition as in (c) using dissipative GRW dynamics. This also leads to anti-thermodynamic behavior. f Comparison of the final temperature distributions of the second (deterministic), third (undamped stochastic), and fourth (damped stochastic) simulation with the initial state of the first one

To test Albert’s proposal, we have performed simulations of a many-particle quantum system with and without GRW collapses. A first aspect to take into account here is the fact that our dynamics should allow for equilibration. This means that it is not possible to use an ideal gas, where “ideal” means that there are no interactions, not even collisions. Such an ideal gas is not guaranteed to approach equilibrium, since all particles will just move on a straight line without affecting each other. As soon as interactions are present, it is not possible to describe the system using a one-particle Wigner function without a closure for the interaction part of the dynamics. However, such closures often already introduce prior assumptions regarding the approach to equilibrium (such as the time-reversal symmetry breaking associated with the hypothesis of molecular chaos, see Sect. 8), which we want to avoid here. Thus, we would have to solve the dynamic Eq. (19) or (20) numerically for a many-particle system over a significant timescale, which is impossible.

We can simplify the problem by noting that we are mainly interested in the effect of GRW collapses on the approach to equilibrium, and not in quantum effects such as superpositions. Using the Wigner framework, we have been able to derive the Fokker-Planck equation (20) which, although describing a quantum system, is very similar in form to classical Fokker-Planck equations. Thus, we can assume that we work in the classical limit where quantum effects apart from collapses (which, in the Fokker-Planck approximation, correspond to momentum diffusion) are negligible, such that particles have well-defined trajectories that can be described using Langevin equations.⁶ The Langevin equations corresponding to the Fokker-Planck equation (20) read

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\vec {r}_i(t)&= \frac{\vec {p}_i}{m}, \end{aligned}$$

(51)

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\vec {p}_i(t)&= - \vec {\nabla }_{\vec {r}_i}U(\{\vec {r}_j\}) + \vec {\chi }_i(t), \end{aligned}$$

(52)

where $\vec {\chi }_i(t)$ is white noise with the properties

$$\begin{aligned} \langle \vec {\chi }_i(t)\rangle _\mathrm {E}&=\vec {0}, \end{aligned}$$

(53)

$$\begin{aligned} \langle \vec {\chi }_i(t)\otimes \vec {\chi }_j(t')\rangle _\mathrm {E}&=2 D_p \mathbbm {1}\delta _{ij}\delta (t-t') \end{aligned}$$

(54)

with the ensemble average $\langle \cdot \rangle _\mathrm {E}$, dyadic product $\otimes$, and identity matrix $\mathbbm {1}$. Thus, we can base our simulations on the Langevin equations (51) and (52). These are numerically much easier to handle than the many-particle Wigner function and thus allow us to model a large system. Moreover, this approach is in the spirit of Albert’s discussion, which is based on the idea that the role of quantum collapses is to provide random perturbations, here represented by the noise term in Eq. (52). If we use a dissipative model such as Eq. (44) instead, Eq. (52) has to be replaced by

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\vec {p}_i(t)= -\gamma \vec {p}_i - \vec {\nabla }_{\vec {r}_i}U(\{\vec {r}_j\}) + \vec {\chi }_i(t) \end{aligned}$$

(55)

and the diffusion constant $D_p$ in Eq. (54) is given by $D_p = \gamma m k_{\mathrm {B}} T_\mathrm {n}$. This then recovers the standard equations of motion for an underdamped Brownian particle, from which the GRW-based model given by Eqs. (51) and (52) can be obtained by taking the limit $\gamma \rightarrow 0$ at fixed $D_p$. We use a smoothed Lennard-Jones potential for the particle interactions (see “Appendix B”). For our simulations, Eqs. (51) and (52) are nondimensionalized (see “Appendix C”) and then solved using the Leapfrog algorithm with time step size $\mathrm {d}t = 0.0025$.

To see how our simulations should be set up, we take into account what precisely Albert’s suggestion is based on: Although, if the initial state of the system is chosen randomly, it is much more likely that the system evolves towards equilibrium, there are certain (“abnormal”) initial conditions for which the system does not evolve towards equilibrium. It is then the (expected) effect of the GRW collapses that the system is brought back onto the right track.

To see whether this actually works, we first require such an abnormal initial condition. For this purpose, we consider a many-particle system whose time evolution is governed by Eq. (8) (i.e., by unitary quantum mechanics without spontaneous collapses). Starting from an initial nonequilibrium distribution at $t=0$, we study its time evolution by solving the unitary dynamics (in its Newtonian approximation given by Eqs. (51) and (52) with $\vec {\chi } = \vec {0}$) numerically. After a certain time $t_\mathrm {rev}$, the system has evolved towards a distribution that is relatively close to an equilibrium state. We then stop and use the distribution at time $t=t_\mathrm {rev}$ with reversed momenta as the initial condition for a second simulation. By the time-reversal invariance of quantum and classical mechanics, the system should then evolve back towards the initial distribution. Since it has thereby moved from a close-to-equilibrium to a far-from-equilibrium state, we have thus found an “abnormal” initial condition (distribution at $t=t_\mathrm {rev}$ with reversed momenta) with the property that the system spontaneously moves away from equilibrium.

Next, we perform another simulation with the same abnormal initial condition, but this time using GRW theory. As is obvious from Eq. (20), the dynamics with spontaneous collapse is no longer time-reversal invariant (the left-hand side and the first two terms on the right-hand side change signs under $(t,p)\rightarrow (-t,-p)$, the third term on the right-hand side does not). If spontaneous collapses can enforce the validity of the second law of thermodynamics in the way suggested by Albert, then the system should be observed to evolve towards equilibrium rather than away from it this time.

The results are shown in Fig. 1, visualizing the local temperature measured by the kinetic energy density $e_{\mathrm {kin}}$ (see “Appendix C”, all variables are dimensionless). We consider a dense fluid (the dimensionless⁷ number density is $\varrho =0.7$), which makes the simulations more efficient and therefore allows us to see all relevant effects. For the first simulation, we prepare two two-dimensional systems of the same size and let them equilibrate at different temperatures $T_1 = 0.5335$ and $T_2=0.6391$.⁸ Afterwards, the systems are brought in thermal contact. This reference state, which has a sharp temperature gradient, is shown in Fig 1a. The total system contains about $8.5\cdot 10^6$ particles. We consider a quadratic domain with boundary length $L=3461.75$.⁹ If we evolve the system forward in time following standard Hamiltonian dynamics (Eqs. (51) and (52) with $\vec {\chi }(t) = \vec {0}$) up to time $t_\mathrm {rev}=625$,¹⁰ we observe—in agreement with what one would expect from thermodynamics—an equilibration towards a state with a more homogeneous temperature distribution. This relaxed state is shown in Fig. 1b. A more detailed analysis of this classical relaxation process can be found in Ref. [52].

For the second simulation, we use as an initial condition the final state of the first one (shown in Fig. 1b), but with reversed particle momenta. This creates an “abnormal” initial condition that leads to anti-thermodynamic behavior under the Hamiltonian time evolution. The simulations confirm this: if we evolve this state forward in time using Hamilton’s equations, we see that the system evolves from a state with a relatively homogeneous temperature (Fig. 1b) towards a state with a sharp temperature gradient shown in Fig. 1c, where the two subsystems have regained their initial temperatures. Consequently, the system has moved away from thermal equilibrium. This is surprising from a thermodynamic perspective, but should be expected from the fact that the underlying Hamiltonian microdynamics is time-reversible. Moreover, effects of this type are known from spin systems [117, 118, 135].

Finally, Albert’s suggestion is tested via the third simulation, which has the same initial condition as the second simulation, but uses the GRW-based model given by Eqs. (51) and (52) with nonvanishing noise $\vec {\chi }(t)$ (see “Appendix C” for the noise amplitude). If his suggestion is correct, then the stochastic perturbations induced by the GRW dynamics should lead to “normal” thermodynamic behavior despite the abnormal initial condition. As can be seen from Fig. 1d, which shows the final state of the third simulation, this is not the case: the system moves away from thermal equilibrium towards a state with inhomogeneous temperature even in the presence of stochastic GRW-type perturbations. We have also performed a fourth simulation with the same initial condition as the second and third simulations, but using the dissipative GRW model given by Eqs. (51) and (55). The final state of this simulation is shown in Fig. 1e, revealing that damping also does not lead to a restoration of thermodynamic behavior. Figure 1f, which compares the temperature profiles for the final states of the second (green), third (orange), and fourth (red) simulation to the initial state of the first simulation (blue), shows that the system evolves back to the initial state (with a small difference due to unavoidable numerical errors) regardless of whether or not stochastic GRW-type perturbations are present. Consequently, our results refute Albert’s proposal according to which, for anti-thermodynamic initial conditions, thermodynamic behavior can be recovered via stochastic GRW-type perturbations.

The numerical experiments presented here are similar to the ones by Orban and Bellemans [136], where anti-thermodynamic behavior was also achieved by time-reversing a molecular dynamics simulation. Orban and Bellemans found that the initial state is not completely recovered due to numerical round-off errors, from which they inferred that anti-thermodynamic behavior is unstable (similar in spirit to Albert’s suggestion). Such “numerical irreversibility” has subsequently been investigated also by other authors (some of whom even linked it to quantum effects) [137, 138]. Since the effects of the round-off errors are now well understood, we can test what happens if GRW noise is added to the dynamics. Figure 1f shows that the difference between the reference state and the result of the deterministic time-reversed simulation (resulting purely from numerical errors) is larger than the difference between the deterministic and the stochastic simulations, even though the amplitude of the GRW noise is much larger than that of the “numerical noise” (such that our “null result” is not a consequence of numerical errors). The temperature gradient in the final state of the time-reversed simulations is smoother than in the reference state, suggesting that the way back from equilibrium has not been completed. However, this increased smoothness is present for both the deterministic and the stochastic simulations. This implies that the GRW noise has no significant effect on whether or not there is anti-thermodynamic behavior, and that the small differences between the initial state of the first simulation and the final states of the second and third simulations are a consequence of round-off errors.

Fig. 2

Time evolution of the absolute values of the Fourier modes of the kinetic energy density $|{\tilde{e}}_{\mathrm {kin}}|$ with wavenumbers a $k=2\pi /L$, b $k=8\pi /L$, and c $k=28\pi /L$ for the first (deterministic, reference), second (deterministic), third (undamped stochastic), and fourth (damped stochastic) simulation. For the second, third, and fourth simulation, time is reversed in this plot to allow for a comparison to the first simulation. The difference between deterministic and stochastic time-reversed dynamics (and between the different forms of stochastic dynamics) is larger for larger wavenumbers, but generally smaller than between the reference and the time-reversed simulations

A similar conclusion can be obtained from Fig. 2, which shows the time dependence of the absolute values of the Fourier modes of the kinetic energy density $|{\tilde{e}}_{\mathrm {kin}}|$ for the wavenumbers (a) $k=2\pi /L$, (b) $k=8\pi /L$, and (c) $k=28\pi /L$ for all four simulations (see “Appendix C”, all variables are dimensionless). In this section and in “Appendix C”, k denotes the wavenumber, and not the temperature parameter as in Sect. 7 and in “Appendix A”. The second, third, and fourth simulations are plotted with reversed time to allow for a comparison to the first one. For all three wavenumbers, the agreement between the two time-reversed simulations is better than their agreement with the first (reference) simulation, implying that the GRW collapses make no significant difference for the macroscopic (thermodynamic) behavior of the system. Therefore, it is unlikely that they are responsible for the emergence of macroscopic irreversibility. However, the collapses do make a certain difference: as can be seen from comparing Fig. 2a–c, the disagreement between deterministic and stochastic simulation results (i.e., the effect of GRW noise) becomes larger for increasing wavenumber k. Since larger wavenumbers correspond to smaller length scales, this means that there are notable differences between deterministic and GRW dynamics on microscopic scales. However, these differences become less important on the macroscopic (thermodynamic) level. A similar result is found when comparing the damped and the undamped GRW model: there are differences on smaller length scales (larger wavenumbers), but not on the macroscopic level.

The results of our simulations thus strongly indicate that, contrary to a widespread claim in the foundations of statistical mechanics, spontaneous wavefunction collapses can not explain the approach to thermodynamic equilibrium even if GRW theory is correct. To understand this result in more detail, we should recapitulate what has led Albert [54,55,56] to the idea that such an explanation would be possible, namely (a) that the property of being an abnormal initial microstate (one that induces anti-thermodynamic behavior) is highly unstable under perturbations, and (b) that GRW theory provides the right kind of perturbations that bring one from an abnormal to a normal microstate. In the simulations, it has been found that the initial macrostate is not completely, but approximately restored, and that the effect of “numerical irreversibility” is much larger here even though the GRW noise has a larger amplitude. This suggests that (a) abnormal microstates are not really stable, but also not as unstable as it is commonly believed, and (b) that GRW theory does not provide the “right” kind of perturbations in this context.

At this point, we should discuss three possible objections to our result. First, it can be argued that the relevant GRW parameters (in particular $D_p$) are so small that observable differences should be expected for macroscopic systems ($10^{23}$ particles), but not for systems with only a few particles (Albert [54] explicitly notes that, if his explanation is correct, an absolutely isolated gas with about $10^5$ particles would not be expected to show thermodynamic behavior). A macroscopic system would be far too large to be studied with microscopic simulations on reasonable timescales. To compensate for this fact, we have used a very strong noise. Assuming, for the sake of definiteness, that our fluid consists of argon atoms, our simulation parameters correspond to $D_p \approx 2.554 \cdot 10^{-43}$ $\hbox {J}^2$s/$\hbox {m}^2$ (see “Appendix C”), whereas using the parameter values $\lambda = 10^{-16}$/s and $\alpha = 10^{14}$ $\hbox {m}^{-2}$ used in Ref. [38] (to get an idea for the order of magnitude of $D_p$) gives (in two dimensions) $D_p = \lambda \alpha \hbar ^2/2 \approx 5.56 \cdot 10^{-71}$ $\hbox {J}^2$s/$\hbox {m}^2$. Adler [43] has suggested values of $\lambda$ that are about 10 orders of magnitude larger, increasing the possible value to $D_p \approx 10^{-61}$ $\hbox {J}^2$s/$\hbox {m}^2$, which is still significantly below our value. If we assume that GRW effects become relevant if $\lambda N$ is of the order 1/s, then for $D_p \approx 10^{-61}$ $\hbox {J}^2$s/$\hbox {m}^2$ (corresponding to $\lambda \approx 10^{-6}$/s) one would require about $10^6$ particles to see GRW effects. Since we have about $10^7$ particles and a noise that is about 9 orders of magnitude larger, there is no reason to assume that our “null result” for effects of GRW collapses on thermodynamics is a consequence of the system size. Current bounds for the GRW parameters are given in Ref. [44].

Second, we have tested a quantum-mechanical theory (GRW) using classical molecular dynamics simulations, motivated by the fact that solving the quantum-mechanical equations for a many-particle system is impossible. This can be justified using the Wigner formalism, which shows that if we ignore all quantum effects except for the spontaneous collapses, GRW theory is equivalent to a classical statistical theory for a system whose dynamics is governed by Langevin equations. Obviously, the result we would have obtained from a full quantum-mechanical calculation would be different. Nevertheless, if a full quantum-mechanical calculation would have led to a different result regarding the problem of irreversibility, then irreversibility would not (solely) be a consequence of GRW collapses, but (also) of the other quantum effects (such as the peculiarities of quantum chaos) that have been ignored in our simulations. Clearly, for a system that is in a macroscopic superposition, the GRW collapses will make an important difference for the dynamics (namely, they will destroy this superposition). However, for explaining why a system of particles that already behave quasi-classically approaches thermodynamic equilibrium, they appear to play no significant role.

Let us consider in particular the position of Wallace [29], who argued that the Wigner function can in principle not be interpreted as a distribution over the positions and momenta of particles in the classical sense, and that even for initially localized particles the classical microdynamics will quickly cease to be accurate in a typical many-particle system due to delocalization and entanglement¹¹. From this view, it might appear inappropriate to simulate a quantum system using classical many-body dynamics. However, one should note here that Albert’s argument (as emphasized above) requires only the fact that GRW theory leads to random perturbations, but not any details of the quantum dynamics. Moreover, the existence of GRW collapses will generally localize the particles and thereby increase the accuracy of the quasi-classical approximation. Finally, even if one denies the existence of localized microscopic particles, solving Eq. (51) and (52) still gives an indication of the solution of Eq. (20) (which is what we are interested in) due to the mathematical relation of Langevin and Fokker-Planck equations.¹²

Third, we have made the approximation of moving from a master to a Fokker-Planck equation [e.g., from Eq. (13) to Eq. (17)] in the derivation of the Langevin equations (51) and (52). This corresponds to the assumption that the Gaussian function that the Wigner function is convoluted with in Eq. (13) decays very rapidly for larger values of $|p'|$. Since the full width at half maximum of this Gaussian is $2\sqrt{\ln (2)\alpha \hbar ^2} \approx 1.756 \cdot 10^{-27}$ kg m/s (which is extremely small), this assumption should also be unproblematic.

11 Conclusion

This article has two central results. First, we have derived master and Fokker-Planck equations for the Wigner function based on the GRW theory, the CSL model, the Diósi-Penrose model, and the dissipative GRW model. This provides a dynamical theory in phase space that allows to explain the emergence of classicality from quantum mechanics regarding both the emergence of Liouville dynamics and the solution of the quantum measurement problem. Moreover, it makes the advantages of the Wigner framework available to researchers interested in spontaneous collapse models. Second, we have used Langevin equations derived from the GRW Fokker-Planck equation for a numerical test of Albert’s suggestion that GRW collapses might be responsible for the approach to thermodynamic equilibrium. The simulations reveal that stochastic GRW-type perturbations do not lead to thermodynamic behavior if it is not already present in the deterministic dynamics. Consequently, our results do not support Albert’s idea.

Appendices

Appendix A: Wigner representation of the dissipative GRW model

Here, we present the derivation of Eqs. (42) and (44) from the main text. For this purpose, we first write the operator (41) in the momentum representation, with momentum coordinates denoted by p and q. We find

$$\begin{aligned} \langle p|{\hat{T}}_\mathrm {d}|q\rangle&= \bigg \langle p\bigg |\frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\int {\hbox {d}p'}e^{\frac{\mathrm {i}}{\hbar }p'{\hat{x}}}e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}\nonumber \\&\quad {\hat{\rho }}\, e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}e^{-\frac{\mathrm {i}}{\hbar }p'{\hat{x}}}\bigg |q\bigg \rangle \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\int {\hbox {d}p'}\nonumber \\&\quad \bigg \langle p - p'\bigg |e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}{\hat{\rho }}\, e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k{\hat{p}})^2}\bigg |q - p'\bigg \rangle \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\int {\hbox {d}p'}\nonumber \\&\quad e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k(p-p'))^2}e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}((1+k)p'+2k(q-p'))^2}\nonumber \\&\quad \big \langle p - p'\big |{\hat{\rho }}\big |q - p'\big \rangle \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\nonumber \\&\quad \int {\hbox {d}p'} e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(4k^2p^2 + 4kp p' - 4 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2 + 4 k p' q - 4 k^2 p' q + 4 k^2 q^2)}\nonumber \\&\quad \big \langle p - p'\big |{\hat{\rho }}\big |q - p'\big \rangle , \end{aligned}$$

(A1)

where $|p\rangle$ and $|q\rangle$ are momentum eigenstates with eigenvalues p and q, respectively. We now use the fact that the Wigner function (2) can also be written as [32]

$$\begin{aligned} W(x,p)= & {} \frac{1}{2\pi \hbar }\int {\hbox {d}p''}\nonumber \\&\quad \bigg \langle {p-\frac{1}{2}p'' \bigg |{\hat{\rho }}\bigg |p+\frac{1}{2}p''\bigg \rangle e^{-\frac{\mathrm {i}}{\hbar }x p''}}. \end{aligned}$$

(A2)

Thus, the Weyl symbol of ${\hat{T}}_\mathrm {d}$ is given by

$$\begin{aligned}&T_\mathrm {d}(W(x,p)) \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{2\pi \hbar }\int {\hbox {d}p'}\nonumber \\&\quad \int {\hbox {d}p''} e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8k^2p^2 + 8 k p p' - 8 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2 + 2 k^2 (p'')^2)}\nonumber \\&\quad e^{-\frac{\mathrm {i}}{\hbar }x p''}\bigg \langle p - p' - \frac{p''}{2}\bigg |{\hat{\rho }}\bigg |p- p'+\frac{p''}{2}\bigg \rangle \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{2\pi \hbar }\frac{1}{\sqrt{2 \pi k^2 r_{\mathrm {c}}^2}} \int {\hbox {d}p'}\int {\hbox {d}p''}\int {\hbox {d}x'} \nonumber \\&\quad \; e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8k^2p^2 + 8 k p p' - 8 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2)}e^{-\frac{(x')^2}{2k^2 r_{\mathrm {c}}^2}}\nonumber \\&\quad e^{-\frac{\mathrm {i}}{\hbar }(x-x') p''}\bigg \langle p - p' - \frac{p''}{2}\bigg |{\hat{\rho }}\bigg |p- p'+\frac{p''}{2}\bigg \rangle \nonumber \\&= \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{\sqrt{2\pi k^2 r_{\mathrm {c}}^2}}\int {\hbox {d}p'}\int {\hbox {d}x'}\nonumber \\&\quad e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8k^2p^2 + 8 k p p' - 8 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2)}\nonumber \\&\quad e^{-\frac{(x')^2}{2k^2 r_{\mathrm {c}}^2}}W(x-x',p-p'). \end{aligned}$$

(A3)

Hence, the dissipative GRW master equation for Wigner functions reads

$$\begin{aligned} \partial _t W(x,p)&= \{H(x,p),W(x,p)\}_\star - \lambda (W(x,p) \nonumber \\&\quad - T_\mathrm {d}(W(x,p))) \end{aligned}$$

(A4)

with $T_\mathrm {d}$ given by Eq. (A3). Equations (A4) and (A3) correspond to Eqs. (42) and (43) from the main text.

Since Eq. (A4) with the Weyl symbol (A3) is extremely complicated, we now derive a simplified form. First, we note that Eq. (A3) can be written as

$$\begin{aligned} \begin{aligned}&T_\mathrm {d}(W(x,p))\\&=\frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{2\pi \hbar } \int {\hbox {d}p'}\int {\hbox {d}p''}\\&\quad e^{-\frac{r_{\mathrm {c}}^2}{2\hbar ^2} (8k^2p^2 + 8 k p p' - 8 k^2 p p' + 2 (p')^2 - 4 k (p')^2 + 2 k^2 (p')^2)} \\&\quad \;\, e^{-\frac{k^2 r_{\mathrm {c}}^2 (p'')^2}{\hbar ^2}}e^{-\frac{\mathrm {i}}{\hbar }xp''}\\&\quad \bigg \langle p - p' - \frac{p''}{2}\bigg |{\hat{\rho }}\bigg |p- p'+\frac{p''}{2}\bigg \rangle . \end{aligned} \end{aligned}$$

(A5)

(Equation (A5) is simply the first line of Eq. (A3).) Motivated by the fact that the standard (nondissipative) GRW model is recovered for $k=0$, we now Taylor expand the integrand in Eq. (A5) to first order in k. This gives

$$\begin{aligned} \begin{aligned}&T_\mathrm {d}(W(x,p))\\&\approx \frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\frac{1}{2\pi \hbar }\int {\hbox {d}p'}\int {\hbox {d}p''}\\&\quad \bigg (1-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8 k p p' - 4 k (p')^2)\bigg )\\&\quad e^{-\frac{r_{\mathrm {c}}^2(p')^2}{\hbar ^2}}e^{-\frac{\mathrm {i}}{\hbar }xp''}\\&\quad \bigg \langle p - p' - \frac{p''}{2}\bigg |{\hat{\rho }}\bigg |p- p'+\frac{p''}{2}\bigg \rangle \\&=\frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\\&\quad \int {\hbox {d}p'}\bigg (1-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8 k p p' - 4 k (p')^2)\bigg )\\&\quad e^{-\frac{r_{\mathrm {c}}^2(p')^2}{\hbar ^2}}W(x,p-p'). \end{aligned} \end{aligned}$$

(A6)

The result (A6) is already reminiscent of Eq. (12) from the standard GRW model. In analogy to the treatment in Sect. 4, we now Taylor expand $W(x,p-p')$ around $p=p'$ in Eq. (A6) and find

$$\begin{aligned} \begin{aligned}&T_\mathrm {d}(W(x,p))\\&=\frac{r_{\mathrm {c}}(1+k)}{\sqrt{\pi }\hbar }\\&\quad \int {\hbox {d}p'}\bigg (1-\frac{r_{\mathrm {c}}^2}{2\hbar ^2}(8 k p p' - 4 k (p')^2)\bigg )\\&\quad e^{-\frac{r_{\mathrm {c}}^2(p')^2}{\hbar ^2}}\\&\quad \Big (W(x,p)-\partial _p W(x,p)p' + \frac{1}{2}\partial _p^2W(x,p)(p')^2\Big )\\&\approx (1+2k)W(x,p) +2 k p\partial _p W(x,p)+ \frac{\hbar ^2}{4r_{\mathrm {c}}^2}\partial _p^2 W(x,p), \end{aligned} \end{aligned}$$

(A7)

where we have dropped terms of order $k^2$ and products of k and $\hbar ^2$ in the last step. After inserting Eq. (A7) into Eq. (A4) and defining $\gamma = 2\lambda k$ and $T_\mathrm {n} = \hbar ^2 /(8 m k k_{\mathrm {B}} r_{\mathrm {c}}^2)$, we finally obtain the Fokker-Planck equation

$$\begin{aligned} \partial _t W(x,p)= & {} \{H(x,p),W(x,p)\}_\star \nonumber \\&\quad + \gamma \partial _p (p W(x,p)) \nonumber \\&\quad + \gamma m k_{\mathrm {B}} T_\mathrm {n} \partial _p^2 W(x,p), \end{aligned}$$

(A8)

which is Eq. (44) from the main text.

Appendix B: Interaction potential

The potential in Eqs. (52) and (55) is (in nondimensionalized form) given by

$$\begin{aligned} U(\{\vec {r}_{k}\}) = \frac{1}{2}\underset{i \ne j}{\sum _{i,j=1}^{N}}U_2(\Vert \vec {r}_i - \vec {r}_j\Vert ) \end{aligned}$$

(B1)

with the pair-interaction potential $U_2$. For $U_2$, we use, following Ref. [139], the smoothed Lennard-Jones potential

$$\begin{aligned} U_2(r) = {\left\{ \begin{array}{ll} 4 \big (\frac{1}{r^{12}} - \frac{1}{r^6}\big ) + C_0 + C_2 r^2 &{} \\ + C_4 r^4 + C_6 r^6 + C_8 r^8 &{}\text { for } r \le r_{\mathrm {c}},\\ 0 &{}\text { otherwise} \end{array}\right. } \end{aligned}$$

(B2)

with $r = \Vert \vec {r}\Vert$, smoothing parameters $C_{0},\dotsc ,C_{8}$, and cutoff length $r_{\mathrm {c}}$. We here use $r_c = 3.5$, implying [139] $C_0 = 7.591 016 534 387 729 7 \cdot 10^{-2}$, $C_2 = -1.858 154 796 630 728 4 \cdot 10^{-2}$, $C_4 = 1.819 594 335 725 356 1 \cdot 10^{-3}$, $C_6 = -8.249 874 769 676 778 6 \cdot 10^{-5}$, and $C_8 = 1.442 811 390 098 910 1 \cdot 10^{-6}$. The form of the interaction potential and the choice of the smoothing parameters ensure that the potential and its first four derivatives vanish at $r = r_{\mathrm {c}}$.

Appendix C: Dimensionless quantities

For the nondimensionalization of our equations, we first write Eqs. (51) and (55) as

$$\begin{aligned} m\frac{\mathrm {d}^2}{\mathrm {d}t^2}\vec {r}_i(t) = -\gamma m \frac{\mathrm {d}}{\mathrm {d}t}\vec {r}_i(t) - \vec {\nabla }_{\vec {r}_i}U(\{\vec {r}_j\}) + A \vec {\xi }_i(t), \end{aligned}$$

(C1)

where A is the noise amplitude and $\vec {\xi }_i(t)$ has zero mean and correlation $\langle \vec {\xi }_i(t)\otimes \vec {\xi }_j(t')\rangle _\mathrm {E} = \mathbbm {1}\delta _{ij}\delta (t-t')$. We then rescale $\vec {r}(t) = \sigma \vec {r}_\mathrm {nd}(t/\tau )$, $t = \tau t_\mathrm {nd}$, $U(\{\vec {r}_j\}) = \epsilon U_\mathrm {nd}(\{\vec {r}_j/\sigma \})$, $\vec {\xi }_i(t) = \vec {\xi }_{i,\mathrm {nd}}(t/\tau )/\sqrt{\tau }$, $\gamma = \gamma _\mathrm {nd}/\tau$, and $A =m\sigma A_\mathrm {nd}/\tau ^{\frac{3}{2}}$, where $\sigma$ is a length scale, $\tau = \sigma \sqrt{\frac{m}{\epsilon }}$ a time scale, $\epsilon$ an energy scale, and the subscript nd denotes a nondimensionalized quantity. (This subscript is omitted in the other sections and the figures.) We then obtain

$$\begin{aligned} \frac{\mathrm {d}^2}{\mathrm {d}t^2}\vec {r}_{i,\mathrm {nd}}(t_{\mathrm {nd}})&=- \gamma _\mathrm {nd}\frac{\mathrm {d}}{\mathrm {d}t}\vec {r}_{i,\mathrm {nd}}(t_{\mathrm {nd}}) - \vec {\nabla }_{\vec {r}_{i,\mathrm {nd}}}U_\mathrm {nd}(\{\vec {r}_{j,\mathrm {nd}}\}) \nonumber \\&\quad \,+ A_\mathrm {nd} \vec {\xi }_{i,\mathrm {nd}}(t_{\mathrm {nd}}). \end{aligned}$$

(C2)

As discussed in Sect. 10, it is useful to compare the strength of the noise used in our simulations to the values expected from current estimates of the GRW model parameters. For the purposes of testing Albert’s theory, it is good if the noise is “too large”, since we can thereby ensure that we are not overlooking any effects due to our “small” system size. An estimate of the order of magnitude is sufficient. The value of $D_p$ in SI units is given by

$$\begin{aligned} D_p= \frac{A_\mathrm {nd}^2\epsilon ^2\tau }{2\sigma ^2}. \end{aligned}$$

(C3)

In the Lennard-Jones potential, the parameters $\sigma$ and $\epsilon$ measure the range and depth of the potential, respectively. To get concrete numbers, we take the example of argon, for which the parameters $\sigma$ and $\epsilon$ take the values $\sigma = 3.4 \cdot 10^{-10}$ m and $\epsilon \approx 1.657 \cdot 10^{-21}$ J, respectively [140]. Moreover, the mass of an argon atom is $m \approx 6.634 \cdot 10^{-26}$ kg [141], such that $\tau \approx 2.151 \cdot 10^{-12}$ s. From Eq. (C3) and $A_\mathrm {nd} = 10^{-4}$ (value used in the simulations¹³), we then find $D_p \approx 2.554 \cdot 10^{-43}$ $\hbox {J}^2$s/$\hbox {m}^2$.

For the simulation with damping, we have used a noise amplitude $A_\mathrm {nd}=10^{-4}$ (as in the undamped case) and a noise temperature $T_\mathrm {n,nd}=0.5863$ related to the noise amplitude via the fluctuation-dissipation relation $A_\mathrm {nd}^2=2\gamma _\mathrm {nd}T_\mathrm {n,nd}$, which determines the dimensionless damping parameter $\gamma _\mathrm {nd}$. For the example of argon, this corresponds to a physical value $T_\mathrm {n} =T_\mathrm {n,nd} \epsilon /k_{\mathrm {B}}=70.356$ K of the noise temperature.¹⁴

The dimensionless temperature $T_\mathrm {nd}$ (we here use $T_1 = 0.5335$ and $T_2=0.6391$) is defined as

$$\begin{aligned} T_\mathrm {nd} = \frac{k_{\mathrm {B}} T}{\epsilon } = \frac{1}{dN}\sum _{i=1}^{N}\vec {v}_{i,\mathrm {nd}}^2 \end{aligned}$$

(C4)

with the dimensional temperature T, number of spatial dimensions d (we use $d=2$ for the simulations), and dimensionless velocity $\vec {v}_{i,\mathrm {nd}}$ of the i-th particle. For argon, the temperatures used here correspond to $T_1\epsilon /k_{\mathrm {B}}= 64.02$ K and $T_2 \epsilon /k_{\mathrm {B}}= 76.692$ K.¹⁵

The local temperature is measured by the local kinetic energy density

$$\begin{aligned} e_{\mathrm {kin}}(\vec {r}_{\mathrm {nd}}) = \frac{1}{2}\sum _{i=1}^{N}\vec {v}_{i,\mathrm {nd}}^2\delta (\vec {r}_{\mathrm {nd}}-\vec {r}_{i,\mathrm {nd}}). \end{aligned}$$

(C5)

Moreover, the y-averaged local kinetic energy density shown in Sect. 1d is given by $\langle e_{\mathrm {kin}}\rangle _y = \frac{1}{L}\smallint {\hbox {d}y_\mathrm {nd}}e_{\mathrm {kin}}(\vec {r}_{\mathrm {nd}})$. For the figures, the local kinetic energy density was calculated via a mass redistribution technique¹⁶.

The Fourier coefficients shown in Sect. 2 are defined as

$$\begin{aligned} {\tilde{e}}_{\mathrm {kin}}(\vec {k})=\frac{1}{2N}\sum _{i=1}^{N}\vec {v}_{i,\mathrm {nd}}^2e^{-\mathrm {i}\vec {k}\cdot \vec {r}_{i,\mathrm {nd}}}, \end{aligned}$$

(C6)

where $\vec {k}= (2\pi /L)\vec {n}$ with $\vec {n}\in {\mathbb {N}}^2$. The definition (C6) ensures that ${\tilde{e}}_{\mathrm {kin}}(\vec {0})=T_\mathrm {nd}$ for $d=2$ [see Eq. (C4)]. In Fig. 2, we always set $k_y = 0$, such that $k = k_x = \Vert \vec {k}\Vert$.