Classical, quantum and event-by-event simulation of a Stern-Gerlach experiment with neutrons

We present a comprehensive simulation study of the Newtonian and quantum model of a Stern-Gerlach experiment with cold neutrons.By solving Newton's equation of motion and the time-dependent Pauli equation, for a wide range of uniform magnetic field strengths, we scrutinize the role of the latter for drawing the conclusion that the magnetic moment of the neutron is quantized. We then demonstrate that a marginal modification of the Newtonian model suffices to construct, without invoking any concept of quantum theory, an event-based subquantum model that eliminates the shortcomings of the classical model and yields results that are in qualitative agreement with experiment and quantum theory. In this event-by-event model, the intrinsic angular momentum can take any value on the sphere, yet, for a sufficiently strong uniform magnetic field, the particle beam splits in two, exactly as in experiment and in concert with quantum theory.


I. INTRODUCTION
In 1922, O. Stern and W. Gerlach demonstrated experimentally that silver atoms passing through an inhomogeneous magnetic field experience deflections in spatially different, distinguishable directions. This observation was very important for the early development of quantum theory for it provided direct experimental evidence that not only the spectra of atoms but also the magnetic moment of the particles might be quantized [1][2][3]. The Stern-Gerlach (SG) experiment is often used in textbooks [4][5][6][7] to introduce the concepts of spin and quantization of angular momentum and plays a prominent role in discussions on determining properties of atomic size objects by means of macroscopic measuring devices [8,9]. The SG experiment, and its conceptually equivalent experiment with single photons passing through a birefringent crystal, are also used in textbooks to illustrate postulates of quantum theory [4][5][6][7].
In short, an SG experiment involves a source of electrically neutral, magnetic particles, collimators, a magnet generating an inhomogeneous field, and a particle detector; see Figure 1 for a sketch of an SG with cold neutrons. Due to the interaction between the magnetic moment of the particle and the inhomogeneous magnetic field, a particle passing through the latter experiences a force that changes the trajectory of the particle. Note that this reasoning is entirely Newtonian, no concept of quantum theory is entering yet.
Assuming (i) uniformly random orientations of the magnetic moments leaving the source and (ii) a sufficiently large uniform magnetic field, the standard classical picture of the magnetic moment as a spinning top leads to the conclusion that there should be no splitting of the beam [5]; however, under certain conditions [11], to be scrutinized in the present paper, the SG magnet splits the particle beam in two, spatially well-separated directions, in agreement with the outcome of the SG experiment. As the amount of deflection is proportional to the magnetic moment, an SG-like apparatus can be used to measure the magnetic moment of nano-size particles [12,13].
As originally conceived, the SG experiment employs electrically neutral particles. Obviously, this begs the question if it would be feasible to perform a similar experiment to observe the spin of say, electrons [14][15][16][17] or ions [18]. Addressing this interesting question is beyond the scope of the present paper, which focuses on the case of electrically neutral particles only.
The deflection in spatially well-separated directions along the direction of the uniform magnetic field is commonly regarded as an experimental proof that the magnetic moment of the particles is quantized [1,2,4,5]. Labeling the distinct beams by a FIG. 1. Diagram (not to scale) of a Stern-Gerlach experiment with cold neutrons, performed by Hamelin et al. [10]. After passing through the collimators, most neutrons travel along the y-direction. The cone indicates the directions in which the neutron may, but not necessarily have to, leave the magnetic field region.
two-valued variable s = ±1/2 and representing the beams by the corresponding state vectors forms the basis for the well-known quantum-theoretical description of the idealized SG experiment [4-7, 19, 20].
The first aim of the present paper is primarily pedagogical in that we present, to the best of our knowledge, the first comprehensive treatment of both the Newtonian and quantum model of a real SG experiment. In order to touch base with a real SG experiment, we have taken model parameters from an SG experiment performed with cold neutrons [10]. In this respect, there is little overlap with earlier numerical studies of the quantum model of an SG experiment [21,22].
The second aim is to demonstrate that a minor modification to the classical, Newtonian equations of motion in the spirit of the event-by-event simulation approach yields results that (i) can be very different from those of the classical and (ii) are in full qualitative agreement with SG experiments and with the quantum-theoretical description thereof. The idea behind this modification is the following. As long as the particle does not experience a magnetic field, the internal frame of reference used to define the direction of the magnetic moment is detached from the laboratory frame of reference. This hold true in quantum theory as well: in the absence of an electromagnetic field there is no relation between the xyz-coordinates of the particle and xyz-components of the spin operator [7]. In the event-based approach, a particle moving from a field-free region into a region where the electromagnetic field is present is viewed as an event which establishes the relation between the xyz-coordinates of the particle and xyz-components of the magnetic moment. This event-triggered process of alignment may be thought of as a highly simplified model for the classical electrodynamic transient processes that occur when a magnetic moment moves through a region in which the magnetic field changes [23].
The paper is structured as follows. Section II describes the SG experiment with neutrons [10] that we take as reference for our simulation work. In Sections III and IV, we present and discuss the results obtained by solving Newton's equation of motion and the time-dependent Pauli equation (TDPE), respectively. Adopting the parameters for the cold neutrons SG experiment in combination with the macroscopic size of the experimental setup requires the use of high-precision solvers and substantial computer resources. In Section IV, we also discuss the transition from a description in terms of position and spin to a model that involves spin-1/2 operators only. Section V introduces the modification to Newton's equation of motion that turns the classical model into a event-by-event, subquantum model for the SG experiment, meaning that data generated by the latter exhibit the same features as the data obtained by SG experiments and their quantum-theoretical description. Section VI summarizes our findings. Figure 1 shows a schematic of the SG experiment with neutrons, as performed by Hamelin et al. [10]. Cold neutrons leaving the neutron guide impinge on a collimator positioned 0.2 m from the exit plane of the neutron guide. The selected neutrons impinge on a second collimator, placed 1 m from the first one. The strongly collimated beam of neutrons then passed through the SG magnet which is 0.8 m long. The distance between the second collimator and the exit plane of the SG magnet is 0.9 m. The direction of the neutrons leaving the SG magnet is selected by means of a meaning window. The distance between the exit plane of the SG magnet and the 3 He detector is 2 m.

II. NEUTRON EXPERIMENT
In Figure 2 we present some of the results reported in Ref. 10. Clearly, the SG magnet causes the neutron beam to split in two well-defined beams, with their maxima of intensities separated by about 6 mm. Note that the window (see Figure 1) in front of the detector moves in the x-direction only.
Looking at the experimental data presented in Figure 2, it is obvious that in order to represent the spin state of a neutron by a two-valued variable, it is necessary to classify the data points as belonging to one of two groups. As the two maxima of the counts are well separated, simply drawing a vertical line at x = 0 suffices to classify the data points. Once this classification is made, we can dispose of the spatial degree of freedom and describe the process of spin-based filtering in terms of spin-1/2 matrices, a model that is often used in textbooks [4][5][6][7].

III. NEWTONIAN MECHANICS
The Hamiltonian describing the dynamics of a neutral, particle of mass m and magnetic moment M subject to a timeindependent magnetic field B = B(x) reads where L is the angular momentum relative to the center of mass x of the particle, and γ is the gyromagnetic ratio. Starting from Equation (1), the standard procedure to derive the equations of motion yields, The angular momentum L has the same dimension ash, namely kg m 2 s −1 . In order to facilitate the comparison with the quantum-theoretical description, it is expedient to define L =hS where S = S(cos φ sin θ , sin φ sin θ , cos θ ) T is a dimensionless vector. In terms of this vector, the classical equations of motion read Note that the presence ofh is the result of rewriting the classical equations of motion in terms of a dimensionless angular momentum S and does not, in any way, imply that Equations (4)-(5) describes quantum phenomena. The length S of the vector S does not affect the solution of Equation (5) and needs to be fixed by comparison with the results of the quantum-theoretical description; this is described in a later section of the paper.

A. Model for the Magnetic Field
Essential to an SG experiment is that the magnetic particles interact with an inhomogeneous magnetic field. Maxwell's equation requires that ∇ ∇ ∇ · B(x) = 0. From the Maxwell equation where J(x,t), ε and µ represent the external current, the electrical permittivity, and magnetic permeability, respectively. It follows that if ∇ ∇ ∇ × B(x) = 0, the magnetic field would induce a nonzero, time-dependent electric field E(x,t). The strength of this electric field would increase linearly with time. Although this electric field would not affect the motion of the electrically neutral particles, in our study, we only consider the case ∇ ∇ ∇ × B(x) = 0. A simple choice, complying with the conditions ∇ ∇ ∇ · B(x) = 0 and ∇ ∇ ∇ × B(x) = 0 just mentioned, is [11,21,22,[24][25][26] that is, B(x) = 0 except when y 0 ≤ y ≤ y 1 where the strength of the field gradient in both the x and z direction is B 1 > 0 (we adopt the convention that B 0 , B 1 ≥ 0). The term in Equation (7) proportional to B 0 , the uniform magnetic field in the zdirection, describes the contribution of the dipole field. The two terms in Equation (7) proportional to B 1 are characteristic for the quadrupole contribution to the magnetic field. The values of B 0 and B 1 depend on the design of the magnet. In this paper, we regard B 0 and B 1 as model parameters.
From Equation (7) it follows that for y ∈ [y 0 , y 1 ], the force F(x) on the particle is given by independent of x or z. For y ∈ [y 0 , y 1 ], the force F(x) on the particle is zero. As a function of y, the simple model Equation (7) shows discontinuities at y = y 0 , y 1 . Instead of smoothing out these discontinuities, we integrate the equations of motion in the interval y 0 ≤ y ≤ y 1 and assume that the velocity distribution at y = y 1 is representative (up to trivial, free-particle scale factors) for the velocity and position distributions at y y 1 . From Equation (8), it follows immediately that the velocity in the y-direction is conserved. In this paper, we assume that all particles move with velocity v y along the y-direction. The time it takes for the particles to traverse the magnetic field region is given by t * = (y 1 − y 0 )/v y .
Once a particle's y-coordinate exceeds y 1 , its velocity v = (v x , v y , v z ) is used to increment the histogram count at the transverse velocity coordinate (v x , v z ) and the simulation of that particle is terminated. The distribution of transverse velocities (v x , v z ) does not change if the particles leave the region where the magnetic field is present and is therefore well-suited to analyze the data. The distribution of transverse positions (x, z) at any plane located to the right of the SG (see Figure 1) is straightforwardly obtained from the distribution of transverse velocities by using the fact that in the field-free region, the particles propagate freely.
In this paper, we mainly present results for the distribution of the transverse velocities (v x , v z ), obtained by classical, quantumtheoretical, and event-by-event simulation. This distribution contains all information about the outcome of the simulated SG experiment and facilitates the presentation of the simulation data in a compact, unified, and convenient manner.

B. Analytically Solvable Cases
It is of interest to consider a special case that is easy to solve analytically. We take as initial positions and velocities of the N particles x = (0, y 0 , 0) and v = (0, v y , 0), respectively, and we only consider the case in which all particles have their initial magnetic moment along the z-axis, i.e., S = S(0, 0, ±1). Note that S × B(x = (0, y, z)) = 0 for any (y, z), see Equation (7), implying that for x = (0, y, z), the torque on the spin is zero; therefore, the direction of the spin does not change and the particles only feel a constant force in the z-direction, see Equation (8). The trajectory is that of a particle in a constant force field, that is v z (t) = ±hγB 1 St/m and z(t) = ±hγB 1 St 2 /2m for 0 ≤ t ≤ t * .
Looking ahead, this simple scenario mimics the quantum-theoretical textbook case (see Appendix D 2) and allows us to fix the magnitude of the classical magnetic moment S. Indeed, the classical and quantum-theoretical expressions for the change in the velocity due to the magnetic field gradients match if S = 1/2. From the analysis of the analytically solvable, classical mechanical case, it follows that the time of flight, the changes of transverse velocity and displacement are given by respectively. The three parameters Equation (9) characterize the state of the particles at the point y = y 1 , that is when they leave the region where the magnetic field is present. Again, looking ahead, the quantum-theoretical textbook case also yields Equation (9). We use Equation (9) to set the scale of time, velocity, and position for both the classical and quantum-theoretical model. The second solvable case is the one that is often referred to when comparing the classical and quantum-theoretical picture of the magnetic moment.
If the uniform magnetic field is present (B 0 > 0), a transformation to a frame rotating with angular frequency γB 0 removes the static field term −γB 0 S z from the transformed Hamiltonian at the cost of introducing time-dependent, sinusoidal terms in the equations of motion. Then, the argument goes, if these sinusoidal terms oscillate sufficiently rapidly, their effect on the motion averages out [6,11]. Although this argument holds for B 0 → ∞, for realistic values of B 0 and B 1 , see Section III C, it does not. Only if the particle trajectories are close to the region where the field gradient is small, the argument applies, see Appendix A. When applied to the SG experiment with realistic values of B 0 and B 1 , the above argument is circular but self-consistent. The justification that the argument is valid comes from the numerical solution presented in Section III E.
If we simply omit the x-component in Equations (7) and (8) (and thereby violate one of Maxwell's equations), we are left with the classical problem in which S z does not change with time and the particle is subject to a force F(x) = γ B 1 S z e z (recall that B 0 has disappeared because of the transformation to the rotating frame). For the initial conditions x = (0, y 0 , 0) and The expressions for final velocity v z (t * ) and position z(t * ) are the same as those obtained in the first analytically solvable case. For each random choice of S, S z is a random number in the range [−1/2, 1/2] and the distribution of velocities is a line at v x = 0, stretching from v z = −v * to v z = v * . This is the expected outcome of the Newtonian description of the SG experiment that is often referred to when comparing with the quantum-theoretical prediction.

C. Model Parameters
We adopt the geometry of the experiment with neutrons, reported in Ref. 10. The region in which there is a nonzero gradient in the x-z directions is 0.8 m long [10], that is y 1 − y 0 = 0.8 m. In the neutron experiment, the maximum gradient of the B-field is estimated to be B 1 = 300 T/m [10]. In the case of the SG experiment with silver atoms, estimates range from B 1 = 1 T/cm = 100 T/m to B 1 = 20 T/cm = 2000 T/m [9,27]. In view of the uncertainties about the strength and precise form of the B-field gradients in these experiments and taking into consideration that the simple form of the B-field gradients used for our theoretical/simulation study is unlikely to hold to any of these experiments, we will use B 1 = 300 T/m in all our simulation work.
In the case of the experiment with neutrons we have [10,28] whereh = 1.05 × 10 −34 kg m 2 s −1 and we have taken as an example B 0 = 1 T. Assume, as we did in Section III B, that all particles have their initial magnetic moment along the z-axis, i.e., S = S(0, 0, ±1). According to Equation (9), the particles cross the plane at y = y 1 at x = (0, y 1 , ±v * t * /2) = (0, y 1 , ±3.53 × 10 −3 m) with a velocity v = (0, v y , ±3.50 m s −1 ). During the remaining free-particle flight to the detector screen, the z-coordinate changes by ∆z screen = ±3.50 × 2 m/395.6 = ±17.7 mm. Thus, in traveling from the source to the detector, the z-coordinate changes by ∆z source−screen ≈ ±21.2 mm. This is about a factor of 7 larger than the splitting observed in the neutron experiment [10]; see Figure 1. In view of the fact that the magnetic field Equation (7) is unlikely to result from the real magnet used in the experiment [10], this order-of-magnitude agreement between the beam-splittings at the screen is quite satisfactory.
As explained in Appendix D, solving the time-dependent Pauli equation for the quantum-theoretical model with the set of parameters given by Equation (10) is computationally very expensive. In order to speed up the development of the simulation software and to generate simulation data for a case that is substantially different than that of neutrons, we have chosen to perform simulations with parameters taken from the original SG experiment [9,27] except that instead of the value of magnetic moment of the silver atom, we have taken the value of the magnetic moment of the Ag 107 nucleus [29]. In the following, we refer to this case as simulations with imaginary silver particles. The parameters are where again, we have taken as an example B 0 = 1 T.

D. Numerical Solution of Equation (III)
In practice, we solve the system Equation (III) by a combination of the exact integration of the torque equation Equation (5) and the velocity-Verlet method as used in molecular dynamics [30]. Appendix B gives the details of the algorithm that we use.
Unless mentioned explicitly, the model parameters for all our classical simulations are B 1 = 300 T/m. Numerical experiments show that the simulation results show insignificant quantitative changes if we decrease the time step from τ = 10 −8 s to τ = 10 −9 s. We use the latter to compute the data that we present in this paper.
Solving Equation (III) for N = 1,000,000 particles with a time step of τ = 10 −9 s takes of the order of hundred minutes on a compute node with two 24-cores Intel Xeon Platinum 8168 CPUs running at 2.7 GHz. We only present data that are essential for the comparison of the classical and quantum description of an SG experiment.

E. Newtonian Dynamics: Simulation Results for Neutrons
In this section, we focus on the SG with neutrons [10]. Repeating the simulations with the particle parameters of imaginary silver particles (see Equation (11) yields data, some of which are presented in Figure 3 and Appendix C, that leads to the same general conclusions.
To allow for a spreading of the particle beam entering the magnet, the distance from the source to the magnet y 0 = 1 m, similar to the distance between the rightmost collimator and the magnet in the neutron experiment [10]. The length of the magnet in the y-direction y 1 − y 0 = 0.8 m [10]. The distance from the magnet to the detection screen is taken to be zero because the motion of the particle is that of a free particle with velocity v = (∆v x , v y , ∆v z ) where ∆v x and ∆v z are the changes of the transverse velocities due to the magnetic field gradients.
As the particles leave the source, the positions and velocities are normally distributed, centered around x = (0, 0, 0) and v = (0, v y , 0) and with variances σ x and σ v , respectively. Unless mentioned explicitly, σ x = σ v = 0.
The simulation reproduces the analytically obtained results if the initial positions and velocities are x = (0, y 0 , 0) and v = (0, v y , 0), respectively, and the initial magnetic moments are aligned along the z-axis, i.e., S = (0, 0, ±1)/2. The results are in excellent agreement with those obtained by solving the problem analytically and are, therefore, not shown.
Next, we assume that the direction of the magnetic moment, represented by the three-dimensional spin vector S, is uniformly distributed over the sphere. In different words, there is maximum uncertainty about the directions of magnetic moments of the neutrons emerging from the neutron guide (see Figure 1). With this initial condition of S, the transverse velocity distribution changes drastically as the strength of the uniform magnetic field decreases from rather strong (B 0 = 1 T) to very weak (B 0 ≈ 0 T), as illustrated in Figure 3a-f.
From Figure 3a, it follows that a uniform magnetic field of B 0 = 1 T is sufficiently strong to suppress the effect of the x-component of the magnetic field. Because the spins S of different particles are distributed uniformly over the sphere of radius S = 1/2, the final distribution of velocities is a strip at v x ≈ 0, stretching from v z = −v * to v z = v * . This is exactly as expected [1,[4][5][6][7]31] on the basis of the arguments discussed in Section III B.
Figures 3b-e demonstrate that the transverse velocity distribution changes drastically each time we reduce B 0 by an order of magnitude. Analytically predicting any of particular shapes shown in Figure 3c-e seems to be a daunting task.
For B 0 = 0, any rotation of the (x, z) coordinates about the y-axis together with the corresponding inverse rotation of the spin leaves the Hamiltonian invariant. As the initial values of (S x , S z ) are distributed uniformly over a circle it follows that the maxima of the transverse velocity distribution are expected to trace out a circle in the v x -v z plane, in agreement with Figure 3f. Figure 4 shows data for the case B 0 = 0. The transverse velocity distribution looks very similar to the one shown in Figure 3f. In the neutron experiment [10], the neutrons that have passed through the SG magnet are selected by means of a narrow window that moves in one direction (say the x direction) only. The recorded neutron counts, plotted as a function of x, show two, very well-separated maxima (see Figure 2). In analogy with the experimental procedure, we compute the one-dimensional, x-dependent distribution by integrating the histogram shown in Figure 4a for v z ∈ [−v * , v * ]/100. This procedure is the computational equivalent of the moving window used in the neutron experiment. The resulting x-dependent distribution is displayed in  In Figure 5, we present the corresponding data for three spin components S x , S y , and S z , obtained by averaging the respective values for v z ∈ [−v * , v * ]/100 (Figure 5a) and v x ∈ [−v * , v * ]/100 (Figure 5b), respectively. Both figures clearly show that the presence of the magnet field gradient causes the initially randomly oriented spins S to preferably align along the direction of transverse propagation. More specifically, focusing on the peaks at v z /v * = ±1 in Figure 5b, we find that the particles with S z ≈ 1/2 (and S x ≈ 0, S y ≈ 0) acquired a negative transverse velocity, whereas those with S z ≈ −1/2 (and S x ≈ 0, S y ≈ 0) acquired a positive transverse velocity, in qualitative agreement with the quantum-theoretical description (see Section IV A). Similarly, looking at Figure 5a, we conclude that particles with S x ≈ 1/2, −1/2 (and S y ≈ 0, S z ≈ 0) acquired a positive (negative) transverse velocity, also in qualitative agreement with the quantum-theoretical description (see Section IV A). The fact that for the x-direction, positive and negative are interchanged with respect to the case of the z-direction is a direct consequence of the different signs of the corresponding components of the magnetic field (see Equation (8)).
Viewed along one direction, e.g., the z-direction, there are two well-separated beams, each of which has a well-defined magnetization. Thus, in the absence of the uniform magnetic field (B 0 = 0), the classical Newtonian model yields a onedimensional profile that displays all signatures of the "quantization of the magnetic moment". Or, put differently, unless the uniform magnetic field B 0 is sufficiently strong, the classical Newtonian model predicts "quantization of the magnetic moment" in any direction.
For completeness, Figure 6 shows how a spread in the initial transverse velocities affects the final transverse velocity distribution for B 0 = 0. Clearly, the main features displayed in Figure 4 are prominently present. Finally, Figure 7 shows that performing the classical simulation using model parameters appropriate for imaginary silver particles instead of neutrons does not change the qualitative features of the transverse velocity distribution. Compared to neutrons (see Figure 4b), the main difference is that the transverse velocity distribution is more spread out over the circle with radius v * (see Figure 7b

IV. QUANTUM-THEORETICAL MODEL
The Hamiltonian describing a neutral, spin-1/2 particle of mass m subject to a time-independent magnetic field B = B(x) reads where p = (p x , p y , p z ) = −ih∇ ∇ ∇ are the momentum operators, σ σ σ = (σ x , σ y , σ z ) are the three Pauli matrices and γ is the gyromagnetic ratio. The magnetic field B(x) is given by Equation (7). From Equations (7) and (12), it follows immediately that [p y , H] = 0, that is, the momentum in the y-direction is conserved. In other words, the motion of the particles in the y-direction is that of a free particle; therefore, in the region where the magnetic field is nonzero, the quantum-theoretical problem effectively amounts to solving the TDPE for the two-component spinor where the subscript s = ±1 refers to the eigenvalues s of the σ z operator.
In Appendix D, we discuss the details of the analytical and numerical tools we use to solve Equation (13).
A. Quantum Theory: Simulation Results Figure 8 shows the transverse velocity distribution | x, z|Φ(t * /10) | 2 obtained by solving the TDPE Equation (D15) for various strengths of the uniform magnetic field and up to the time t * /10 at which, in the Newtonian model, the neutrons would have left the region in which the magnetic field is present.
For a sufficiently strong uniform magnetic field, e.g., B 0 = 1 T, the transverse velocity distribution is bimodal with wellseparated maxima at v z ≈ ±v 0 ; see Figure 8a. The SG magnet then functions as an (almost perfect) filtering device, yielding particle beams which may be labeled by the eigenvalues of the σ z Pauli matrix.
We wrote may because a meaningful assignment in terms of the eigenvalues σ z requires that if we send the beam of particles through a second SG magnet with its strong uniform magnetic field along the z-axis, the particles should emerge in one and the same beam only.
More generally, if we use a filter device to label different outcomes, subsequent repeated filtering by identical devices should leave the labeling intact [32]. If it does not, the original assignment is useless.
Thus, to verify that an SG magnet with its strong uniform magnetic field along the z-axis acts as a spin-filtering device, we repeat the simulation with B 0 = 1 T and initial spin state | ↑ . The resulting transverse velocity distribution is the same as the one in Figure 8a with the top spot removed (image not shown). Thus, with a strong static field B 0 , the SG magnet indeed acts as an ideal filtering device.
In the quantum-theoretical treatment, the spin is quantized by construction; therefore, the observed splitting of the beam cannot be regarded as evidence for the quantization of the spin; however, for large B 0 , the quantized spin model shows that the SG magnet splits the beam (in agreement with experiment) whereas the Newtonian model does not (in disagreement with experiment), exposing a fundamental shortcoming of the latter.
As in the classical case (see Figure 3a- Maxwell's equation dictates that (with our choice of the frame of reference) the Hamiltonian of an SG experiment should contain terms in both γσ x B 1 and γσ z B 1 , which implies that the magnetization (in any direction) is not conserved. Therefore, unless B 0 → ∞, the eigenvalues of σ z cannot be used to label the eigenstates of the Hamiltonian. In other words, there are situations, choices of the model parameters, for which the SG magnet cannot be used to define the quantization direction of the spin [21,22].
We study this aspect by solving the TDPE for the initial state given by Equation (D18) with θ = α = 0, that is for the initial spin states | ↑ and | ↓ and B 0 = 0. The transverse velocity distributions are shown in Figure 9a,b. If the initial spin state is | ↑ (| ↓ ), the wave packet dominantly propagates along the −z-direction and +z-direction; see Figure 9a,b, respectively. The term Not surprisingly, the sum of Figure 9a,b yields an image that looks very much like Figure 8f. Figure 10a,b shows the corresponding probability distributions for the | ↑ and | ↓ components of the wave function, projected onto the z = 0 and x = 0 axis, respectively.  Figure 9a. The probability distributions |Φ −1 (v x = 0, v z ,t * /10)| 2 is too small to be visible in the plot. Except for |v x /v 0 | ≤ 0.2, the difference between two distributions is too small to be visible in the plot. For presentation purposes, each distribution is normalized such that its maximum is one. As in Figure 9, v 0 = v * /10.
If, in an experiment such as the one with cold neutrons [10], one would only count particles by moving a narrow window along the x-direction, the distribution shown in Figure 10a would lead us to conclude that the SG magnet has split the beam into parts. On the other hand, measuring with a moving window along the z-direction yields the distribution shown in Figure 10b, which forces us to conclude that only the | ↓ component is present in the outgoing beam. Indeed, the intensity of the | ↑ component is several orders of magnitude smaller than the one of the | ↓ component. For B 0 > ∼ 0, the SG magnet does not act as a spin filter.
It may be of interest to note that if an SG magnet is used to measure the magnetic moment to, e.g., atomic clusters [13], the value of B 0 does not matter much. The positions of the peaks in the one-dimensional distributions, which are the same for large and zero uniform magnetic field B 0 , suffice to determine the value of the magnetic moment.

B. Quantum Theory: Simplified Model
The TDPE Equation (13), with the term in σ x removed, is an excellent approximation to the full TDPE Equation (13) if the uniform magnetic field is strong enough, e.g., if B 0 = 1 T. Then, we may also replace σ z by σ σ σ · b where b = B/ B is the unit vector parallel to the strong, uniform magnetic field, because (i) the eigenvalues of σ σ σ · b are the same as those of σ z and (ii) only the eigenvalues enter in the spin part of the simplified TDPE. By introducing b, the latter can describe situations in which the strong uniform magnetic field can take any orientation, as long as it is approximately perpendicular to the y-direction (otherwise the argument to remove the σ x term may break down).
We can now simplify the description further. Because of the one-to-one correspondence between the eigenvalue of σ σ σ · b and the change in the transverse velocity of the outgoing particles, we may dispose of the description of the translational degrees of freedom entirely and represent the operation of the SG apparatus by the projection operator [7] acting on the spin state |ψ = a ↑ | ↑ + a ↓ | ↓ only. The probability to observe a particle in the beam labeled by s b = ±1, one of the two eigenvalues of b · σ σ σ , is given by where cos ξ = s · b and s = ψ|σ σ σ |ψ . The last expression in Equation (16) is reminiscent of Malus' law for the intensity of polarized light passing through a polarizer. The projector equation (Equation (15)) and the probability equation (Equation (16)) describe the operation of the SG apparatus in terms of the spin-degree of freedom only. This simplified model is often used in textbooks to elucidate quantum measurement theory [5][6][7]. We stress that Equation (16) does not apply to the case of a weak uniform magnetic field.
Solving the TDPE Equation (13) for B 0 = 1 T and for the initial states Equation (D18) with θ = 0, π/6, π/4, π/3 and α = θ /2 yields the expected bimodal shape of the transverse velocity distributions (data not shown). The total probabilities for v z < 0 and v z > 0 are in excellent agreement with the prediction based on Equation (16). The qualitative difference between the Newtonian and quantum-theoretical prediction in the case of a large uniform magnetic field has been decisive to eliminate the former as a description of the experimental observations [1][2][3]; however, that does not imply that quantum theory is the only viable description of experiments in which the frequency distribution of detection events is built up one-by-one, such as in the SG experiment.
From this broader perspective, the fundamental question to be answered is "is it possible to construct a process that generates event-by-event and without using knowledge about the final distribution of events, frequency distributions that are commonly thought to be a signature of wave interference, two-particle entanglement, uncertainty, etc." This question is answered in the affirmative by the event-by-event simulation approach developed in Refs. 33-44. In the case at hand, the conceptually interesting question is whether it is possible to retain a picture of the SG experiment in which individual particles follow trajectories while, in contrast to the Newtonian results, the transverse velocity distribution exhibits two well-separated maxima along the line defined by the direction of the strong static field (the z direction in our case).
Remarkably, a marginal modification of Newton's equation of motion suffices to answer this question affirmatively. The modification consists of replacing step in which the force F is being calculated (see Appendix B for details) by the rules 4. If y ∈ [y 0 , y 1 ]: the first time that the event B(x) > 0 occurs, that is when the particle enters the region where B(x) > 0, use Equation (16) with s = S to align the vector S along the magnetic field s b B(x) and compute F = γ B 1 S z e z − γ B 1 S x e x .
In detail, if r ≤ S · B(x)/ B set S = B/2 B , otherwise set S = −B/2 B . Here r is a uniform (pseudo) random number in the range [−1/2, 1/2] (which changes each time before it is used). With the new S, compute F = γ B 1 S z e z − γ B 1 S x e x . For each particle, the alignment of S is carried out only once.
One might try to argue that because the event-by-event model makes use of Equation (16), it implicitly "knows" about quantum theory; however, probabilistic laws such as Equation (16) also follow from the application of logical-inference [19,45] to the modeling of event-based processes. This approach yields Equation (16) directly, without any reference to quantum-theoretical concepts.
In short, the key idea of the logical inference approach is that "good" physics experiments must yield reproducible frequency distributions which are robust, meaning do not change much, if the conditions under which the data was taken changes a little [45]. In the case at hand, the frequency distribution consists of the average numbers of +1 and −1 events and ξ = arccos(2S · B(x)/ B ) represents the condition [45]. Expressing the key idea mathematically leads to the requirement that the Fisher information for the probability p(x|ξ ) to observe the event x = ±1 under the condition ξ must be independent of ξ and minimal [45]. After some elementary algebra, we find that the solution of this optimization problem reads [45] where the ± sign reflects the ambiguity in assigning +1 or −1 to one of the directions. Quantum theory postulates Equation (16) (through the Born rule) whereas the logical inference approach allows us to derive Equation (16) without making reference to a concept of quantum theory. Thus, the argument that the event-by-event algorithm implicitly refers to quantum theory does not hold. Moreover, the modification does not change the vector character of S. In the event-by-event model, S can take any value on the sphere of radius 1/2, there is no wave function, there are no Pauli spin matrices, there simply is no element of quantum theory in the event-by-event model. Figure 11a demonstrates that the event-based model produces a bimodal transverse velocity distribution, in qualitative agreement with the solution of the TDPE Equation (13). Clearly, the minor modification to Newton's equation has a tremendous impact on the trajectories of the particles. For B 0 ≈ 0, the event-by-event simulations yields the circular distribution; see Fig-ure 11e,f, in qualitative agreement with both the Newtonian and quantum-theoretical description.

VI. CONCLUSIONS
In all our simulations, the strength of the magnetic field gradient was fixed and tuned to the case of an SG experiment with cold neutrons [10] while the strength of the uniform component of the magnetic field was varied. The simulation data for imaginary silver particles instead of neutrons show the same qualitative features.
In Table I, we collect the most essential features of the results for the transverse velocity distribution obtained from the simulation of three different models of the SG experiment. Thereby, we have omitted many of the computational details mentioned earlier and limit the discussion to the two extreme cases of a strong and zero uniform magnetic field.  Fig. 3(f)) one stripe (Fig. 3(a)) Quantum theory circular ( Fig. 8(f)) two spots ( Fig. 8(a)) Event-by-event circular ( Fig. 11(f)) two spots ( Fig. 11(a)) The first three rows of the last column in Table I express what is known since the original SG experiment was performed, namely that Newtonian mechanics cannot explain the observed splitting of the particle beam [1,[4][5][6][7]31] if the uniform magnetic field component is sufficiently large. It is exactly under this last condition that the quantum-theoretical textbook model provides an accurate description of the time evolution of the probability distribution while the Newtonian model does not. However, we have also shown that a minor modification to Newton's equations of motion yields results that are in line with the experimental observation and quantum theory. In this event-by-event simulation approach, the spin is described in terms of a three-dimensional vector, not in terms of Pauli matrices.
If the strength of the uniform magnetic field B 0 gradually decreases, then, for any of the three models, the changes in the transverse velocity distribution become hard to predict analytically, unless the effect of B 0 becomes negligible. Indeed, for B 0 = 0 a symmetry argument can be used to understand why the calculated transverse velocity distribution shows a circular structure, see column two of Table I. However, also the case B 0 > ∼ 0 poses some interesting interpretational issues, depending on how the distribution of particles is measured. If, as in the neutron experiment [10], one only records the distribution along a particular direction, the Newtonian model also yields a bimodal distribution. Without additional data, the bimodality would (erroneously, see Section IV A) imply that the two beams can be labeled by the spin quantum number.
From a general perspective, quantization (not to be confused with results from quantum theory) is the process of classifying empirical data into groups and attaching discrete labels to these groups. As mentioned at the end of Section II, in the specific case of the neutron experiment it is clear that quantization is the result of classification, putting data points in two groups; see Figure 2. Once this "operation" has been carried out, the compressed, new data are "quantized". In our view, quantum theory provides a powerful mathematical framework to describe such "quantized data". Within quantum theory, the spin is quantized by definition/construction. If the "quantized" form of the empirical data is described well in terms of a quantum spin model, then that is a great achievement; however, this success does not necessarily justify the conclusion that "quantization" is a property/attribute of the phenomenon that gave rise to the empirical data. In our view, drawing this conclusion mixes up the phenomenon that gave rise to the empirical data with a quantum model of it.
On the basis of SG experiments that have been performed to date, it is not possible to distinguish between the quantumtheoretical and event-by-event model. New, high-precision experiments are needed to rule out the latter and to allow for a quantitative comparison between experimental and simulation data.
Furthermore, it would be of interest to perform an SG experiment in which the uniform magnetic field is weak enough to render the description textbook model invalid. For instance, an experiment with neutrons passing through a quadrupole magnet with a large field gradient would allow a direct comparison with our simulations (which, if needed, can easily be adapted to other field configurations). and integrating over time gives where β =hB 1 /mB 0 . Integrating Equation (A4) over time once more results in where β = β /γB 0 . In the case of neutrons and for B 0 = 1 T and B 1 = 300 T/m we have β ≈ 2 × 10 −5 m/s and β ≈ 10 −13 m and it follows immediately from Equation (A5) that the motion of the spin has a negligible effect on the motion of the particles in the xdirection.
On the other hand, under the same conditions, it follows from Equations (A1c) and (A2) that S z (t) ≈ S z (0) and the equation of motion for the z-component of the velocity becomes yielding In our simulations, S z (0) is a uniform random number in the range [−1/2, +1/2]. Therefore, if v z (0) = 0, the values of v z (t * ) are also uniformly distributed over the interval [−v * , +v * ].
Clearly, these elementary calculations yield results which are in excellent agreement with the simulation data for B 0 = 1 T.

Appendix B: Numerical Solution of Equation (III)
For any time step τ, Equation (5) can be solved in closed form. In terms of the three components of the spin vector S, we have where where u = γB x /Ω, v = γB y /Ω and w = γB z /Ω, C = cos(τΩ), S = sin(τΩ) and Ω = |γ|(B 2 x + B 2 y + B 2 z ) 1/2 . The matrix R(τ) is orthogonal, implying that the integration scheme does not change the length of S.
We integrate the equations of motion Equation (III) using the velocity-Verlet algorithm [30]. Initially, the positions and velocities are normally distributed, centered around x = (0, 0, 0) and v = (0, v y , 0) and with variances σ x and σ v , respectively.
We only consider the case 0 < y 0 < y 1 and y 0 ≤ y ≤ y 1 . According to Equation (8), at t = 0, the force F(x,t = 0) = 0. We choose a time step τ and repeat steps 1 to 5: for a number of time steps N.
As a curiosity, it may be of interest to mention that if the force F is constant, the Verlet scheme integrates the equation of motions exactly. On the other hand, Equations (B1) and (B2) integrate the torque equation for spin exactly. Thus, it is only the combination of particle and spin motion that forces us to integrate Equation (III) numerically. Appendix C: Newtonian Dynamics: Imaginary Silver Particles Figure 12c demonstrates that a uniform magnetic field of B 0 = 0.01 T is sufficiently large to yield a distribution of velocities centered at v x = 0 and stretching from v z = −v * to v z = v * , as expected for the classical, textbook model [1,[4][5][6][7].
The transverse velocity distribution changes drastically if B 0 decreases two and three orders of magnitude; see Figure 3d,e. As B 0 vanishes, the transverse velocity distribution becomes a circular disk of radius v * , with the maximum located at the edges and the minimum at the center, qualitatively similar to the case of neutrons.
where s = ±1 denotes the eigenvalues of σ z and we omitted the constant term proportional to k 2 x . We eliminate the term proportional to B 0 by a transformation to a rotating frame. Substituting ϕ s (k z ,t) = e istγ B 0 /2 φ s (k z ,t), Equation (D5) becomes Changing to new variables defined by k z = k + gt and t = t , we have ∂ ∂ k = ∂ ∂ k z , ∂ ∂t = ∂ ∂t + g ∂ ∂ k z , and Equation (D6) changes to Choosing the solution of which reads or, in terms of the original coordinates, It then follows that for any choice of the initial state ϕ s (k z , 0). Equation (D11) tells us that as a function of time, the probability density |ϕ s (k z ,t)| 2 is the same as the initial probability density, translated in momentum space by −sγ B 1 t/2. Assuming that ϕ s (k z , 0)|ϕ s (k z , 0) = 1, it follows from Equation (D10) that k z (t) s = ϕ s (k z ,t)|k z |ϕ s (k z ,t) = ϕ s (k z ,t)|(k z − gt)|ϕ s (k,t) + gt ϕ s (k,t)|ϕ s (k,t) = ϕ s (k z , 0)|k z |ϕ s (k z , 0) + gt = k z (0) s + gt = k z (0) s + sγ B 1 t 2 .
(D12) Therefore, in the textbook case, the presence of a gradient in the magnetic field causes the average momentum to linearly decrease (s = +1) or increase (s = −1) if γ < 0 (which is the case for neutrons or imaginary silver particles). Integrating Equation (D12) with respect to time yields z(t) s = q z (t) s = z(0) s + k z (0) s t + sγ B 1 t 2 /4 showing that the average position traces out a parabolic trajectory [7]. Writing Equation (D12) in terms of the quantum-theoretical velocity operator defined by v =hk/m, we have v z (t) ±1 = v z (0) ±1 ±hγ B 1 t/2m. The classical mechanical expression ±hγB 1 St/m for the change of the velocity due to the magnetic field gradient matches the quantum-theoretical result if S = 1/2, as mentioned earlier.

Dimensionless Form
We definehk x /m = v 0 v x ,hk z /m = v 0 v z , and t = t 0 τ where v 0 and t 0 set the scale of the velocity and time, respectively. In terms of these variables Equation (D3) reads We simplify the notation somewhat by introducing the (dimensionless) parameters and, at the risk of creating confusion, make the replacements τ → t, v x → x, v z → z, i ∂ ∂ v x → p x and i ∂ ∂ v z → p z . Then, Equation (D13) becomes i ∂ ∂t |Φ(t) = a x 2 + z 2 − cσ z − bσ z p z + bσ x p x |Φ(t) .
We emphasize that from now on, whenever we discuss the quantum-theoretical model, x and z denote the dimensionless velocity in the xand z-direction, respectively. In the case of neutrons, we use the parameters given in Equation (10) and the corresponding values of t 0 = t and v 0 = v to find whereas for c (using Equation (11)), we find Comparing Eqs. (D16) and (D17), we may expect that numerically solving the TDPE Equation (D15) for the case of neutrons is much more difficult than for the case of imaginary silver particles, simply because in the former a and b differ by more than five orders of magnitude.

Simulation Method
In practice, we solve Equation (D15) by means of the Chebyshev polynomial algorithm [47,48]. Disregarding the discretization of the p x and p z , this algorithm has the virtue that it yields numerical results with close to machine precision, for any time t [47,48]. Technical details about this TDPE solver can be found in Appendix D 6.
For the reasons explained in Appendix D 6, simulating the neutron experiment [10] is prohibitively costly. A simple way out of this conundrum is to "make the SG magnet shorter" by a factor of ten. Then t 0 = t /10, v 0 = v /10 and, according to Equation (D14), a is reduced by a factor of thousand and b is left unchanged.
A simulation with t 0 = t * /10 and v 0 = v * /10, using a grid of 2 15 × 2 15 points requires about 200 GiB of memory and finishes in about 8 hours, using a compute node with two AMD EPYC 7742, 2 × 64 cores, 2.25 GHz processors. Although it would be technically straightforward to use more nodes by extending the code to use the MPI communication protocol, we believe that this extension would not bring much additional insight. The main point here is that we can perform simulations in which the neutrons travel a macroscopic distance and the dimensions of the apparatus are also macroscopic.
It may be of interest to mention that in practice, a product-formula approach [22,49] using the decomposition in terms of the exact propagators in Equation (D1) for the x and z components fails, simply because of the large disparity between the coefficients a and b (in the case of neutrons).