Quantum Dynamical Simulation of a Transversal Stern--Gerlach Interferometer

Originally conceived as a gedankenexperiment, an apparatus consisting of two Stern--Gerlach apparatuses joined in an inverted manner touched on the fundamental question of the reversibility of evolution in quantum mechanics. Theoretical analysis showed that uniting the two partial beams requires an extreme level of experimental control, making the proposal in its original form unrealizable in practice. In this work we revisit the above question in a numerical study concerning the possibility of partial-beam recombination in a spin-coherent manner. Using the Suzuki--Trotter numerical method of wave propagation and a configurable, approximation-free magnetic field, a simulation of a transversal Stern--Gerlach interferometer under ideal conditions is performed. The result confirms what has long been hinted at by theoretical analyses: the transversal Stern--Gerlach interferometer quantum dynamics is fundamentally irreversible even when perfect control of the associated magnetic fields and beams is assumed.


Introduction
The experiment by Stern and Gerlach performed in 1921 (a centennial this year) [1,2] proved to be of fundamental importance for the development of quantum mechanics. Far from being a thing of the past, the original design of the Stern-Gerlach apparatus (SGA) has continued to inspire new questions [3][4][5][6][7][8][9][10][11][12], and experimental techniques [13][14][15][16][17]. Professor David Bohm in his well-known textbook on quantum theory discusses a device consisting of two SGAs (see fig. 1) joined in an inverted manner so that If the uniform magnetic fields (. . . ) are set up in exactly the right way, and if the second inhomogeneous field is an exact duplicate of the first one, the two wave packets can be brought together into a single coherent packet. Although the precision required to achieve this result would be fantastic, it is, in principle, attainable.
-Section 22.11 in [18] 1 The analysis of this possibility gave rise to what is known as the Humpty-Dumpty problem [20][21][22] of coherent recombination of spatially separated partial atomic beams (i.e. the spin-up and spindown components of the spinor wave function) coming out of an SGA, as well as to a plethora of experimental techniques within the domain of matter-wave interferometry, all of which use different approaches to work around the fantastic precision required to realize the original 1 In Bohm's proposal the atom is deflected in regions of "uniform magnetic fields" -but there are no such Lorentz-type forces acting on electrically neutral atoms. All subsequent works, including ours, consider magnetic gradients to achieve the restoring deflection, with three more SGAs completing the interferometer sketched in fig. 1. There is also Wigner's scheme, in which the inhomogeneous magnetic field from an electric current recombines the beams that emerge from the SGA [19]; it has not been used for a quantitative model. . Four Stern-Gerlach magnets are used to generate a magnetic field B that first splits and then recombines the atomic beam of silver atoms, which are moving with velocity v y through the apparatus. The spin quantization axis is determined by the homogenous field B 0 oriented along the Z axis.
In this work we revisit the question of time-reversibility of a transversal SGI. Past theoretical work focused on the impossibility of perfect control of magnetic fields and beams, indicating a quick loss of spin coherence if such control could not be maintained [20,21]. Yet, what was also suggested by these analyses is that the very dynamics of a beam-splitting Stern-Gerlach apparatus might not be reversible as imagined, even when perfect control over involved magnetic fields and beams is assumed [27]. To fill in the missing piece in the quantum Humpty-Dumpty riddle, we have performed an accurate three-dimensional wave-propagation simulation of a transversal SGI in order to capture accurately all quantum mechanical effects that arise once a full, ideal magnetic field is accounted for. The results of our simulation confirm the intuition that the quantum dynamics of a Stern-Gerlach apparatus is not reversed by the application of inverse magnetic fields, even if one has perfect control over the experimental apparatus.

Modeling the Transversal SGI
A beam of silver atoms prepared in a pure spin state directed along the X axis, |ψ = ψ(x, y, z)(|↑ + |↓ )/ √ 2, with the initial wave function ψ(x, y, z) chosen to be a gaussian wave packet of width δ, The atoms enter the apparatus from the left at time t = 0, at x = 0, y = y 0 = −L/2, z = 0, travel in the Y direction with velocity v y , and exit at y = L/2, with y = ±L/2 sufficiently far from the SGI center at x = y = z = 0 to be outside the magnetic fringing fields. Inside the apparatus, the Magnetic field model used in our simulation. In the Y,Z plane we have a fictitious magnetic line charge λ(y) at z = −a and the opposite charge −λ(y) at z = a. For points on the Y axis, the resulting magnetic field is in the Z direction, B = b 0 (y) e z . The model is defined by specifying b 0 (y); see section 2.2. magnetic field B is generated by two symmetrically placed linear "magnetic charge" 2 distributions of opposite charge at the height a from the beam's initial position, according to the arrangement depicted in fig. 2. The resulting magnetic field is symmetric across the XY plane and the nonlinear dependence on the z coordinate of the magnetic field due to a single charge line, which ordinarily would result in slightly different trajectories for the spin-up and spin-down components, no longer poses an obstacle for spatial beam recombination. Furthermore, to fix the spin quantization axis along the Z axis, a bias field B 0 = B 0 e z in the Z direction (here e z is the unit vector in the Z direction, and B 0 is the magnitude of the bias field) is applied throughout the apparatus. The magnitude of this field is chosen to be much greater than the magnitudes of the B x , B y components in the apparatus to suppress the effects of Larmor precession around axes other than Z.
The Larmor spin precession angle that the particle accumulates as it travels through the apparatus results in the final spin of the atom not necessarily pointing in the X direction despite our best efforts, even if the beams are recombined coherently. It is therefore inappropriate to check whether the final spin state is equal to the initial one. Instead we measure the degree of preserved coherence, C(t), by the purity of the statistical operator M(t) for the spin degree of freedom, where σ is the vector of Pauli spin matrices. The spin coherence C(t) ranges between 0.0 and 1.0, between a completely mixed spin state (coherence lost) and a pure spin state (coherence preserved), with intermediate values indicating some degree of coherence loss. After time T, when the particle has traversed the length of the apparatus L and reached the position y = L/2, it emerges at the right and its spin coherence is measured. It is our goal that the final value C(T) of the spin coherence is as close to its maximal value 1.0 as possible; in a perfect experiment we shall have C(T) = 1.0. During the evolution the peak separation between the partial beams reaches ∆Z(T/2) (henceforth denoted by ∆Z). One has to ensure that ∆Z is large compared with the width δ of the atom wave packet to say that the partial beams have truly separated. If despite partial beam separation in the middle of the apparatus the spin coherence has been preserved, then the evolution in the left half of the apparatus has been reversed, or in the playful language of [20], Humpty-Dumpty, which is really an egg, has been put back again after his great fall. After this overview of the experiment we move on to describing the methods used in the simulation.

Quantum Dynamics of the Apparatus
Following the usual quantum treatment of an SGA, we have the following Hamiltonian describing the spin-1/2 silver atom comprising the beam: with the momentum operator P = (P x , P y , P z ), the position operator R = (R x , R y , R z ), the magnetic moment µ, the vector of Pauli spin matrices σ = (σ x , σ y , σ z ), and the position-dependent magnetic field B( R) in the apparatus. The numerical values for m and µ that are used in the simulation, and other simulation parameters are included in the appendix.
To solve the equations of motion with an accuracy high enough to model coherent beam recombination we employ a fourth-order Suzuki-Trotter operator splitting method [30]. In this scheme, given any Hamiltonian with a potential-energy term V( R), the time propagator , which advances the wave function by a time step t, is split into several exponential terms, each involving only the momentum operator or the potential operator with coefficients specially chosen, such that the whole formula is accurate to a certain order k, where A i , B i correspond respectively to potential and momentum operator terms. While approximations U N for increasing N yield higher-order formulas and can in principle achieve any desired level of accuracy, the computational cost goes up considerably with increasing N, as between each application of operators either in momentum or in position space one has to perform a Fourier transformation on the wave function. For our purpose we have chosen a fourth-order Suzuki-Trotter method of a special form that involves a gradient term of the potential and five operator terms in total, compared with 11 terms in the standard Suzuki-Trotter formula needed to achieve the same order [31][32][33]. The time propagator in this "7 approximation" is 3 and to propagate the wave function the above five unitary operators are applied in their respective space. For instance, the application of the first two terms from the right on ψ( R, t 0 ) as part of the algorithm requires two Fourier transformations F , A single application of the U 7 propagator requires the execution of four Fourier transformations for a spinless wave function. Time propagation of a spin-1/2 particle will require eight Fourier transformations per time step, an operation that contributes heavily to the computational cost. Once this split-operator machinery is in place we can simulate the quantum evolution of any system described by a reasonable Hamiltonian with a position-dependent potential, including the quantum dynamics of an SGI. However, when attempting to do so according to the method presented above, one is quickly faced with rising computational space and time resources required to perform the computation to a sufficient degree of accuracy. The storage cost of a wave function across the region of the apparatus is prohibitive, as is the time required to perform repeated Fourier transformations in each time step.
To make the calculation feasible, we transform the Hamiltonian into the co-moving frame of the atomic beam traveling with velocity v y along the Y axis of the apparatus. This is done at the cost of making the potential term in the Hamiltonian time-dependent, and thus requires a re-evaluation of the split-operator method which, presented in its most basic form, is applicable only to Hamiltonians not explicitly dependent on time. The extension of the split-operator method to time-dependent non-commuting Hamiltonians is treated in full generality in [35]. Here we just note the elegant result: in the time-dependent split-operator method U 7 T the reference time has to be advanced cumulatively at each encountered exp − itβ h P 2 2m term, so that subsequently applied potential operator terms use the newly advanced time T , until the next application of a position propagator and so on. The resulting formula for the time-dependent propagator is then guaranteed to be accurate to the same order as the original formula. Using this simple prescription we readily obtain the time-dependent propagator U 7 T , for the evolution from time T to time T + t, In this way one has to keep track only of the wave function in the region required to model the beam deflection in the Z axis, which is very small compared with the total extent of the apparatus, as the movement along the Y direction has been taken care of by transforming the Hamiltonian into the co-moving frame. The outstanding challenge is an efficient and accurate computation of the magnetic field B(x, y, z), which in each time step has to be evaluated across the domain of the wave function, now in the co-moving frame of the beam.

Magnetic Field in the Apparatus
To our best knowledge all previous quantum treatments of Stern-Gerlach-like experiments employed a magnetic field that either did not satisfy the Maxwell equations through the neglect of the X and Y field components and linear approximations [20,21,36], or used a simple model of a magnetic field that neglects the fringing fields at the entrance and exit of the apparatus [37,38]. For the purpose of our work it is crucial that the employed field model is free from both of these approximations, since it is clear that either of these simplifications will result in an apparatus that will seemingly have no problem with reuniting the partial beams in a coherent manner, as long as the beams and fields are accurately controlled.
In the traditional design, at least four Stern-Gerlach magnets have to be employed. As the particle goes through the apparatus, the partial beams corresponding to the spin-up and spin-down components are first imparted a momentum which leads to their spatial separations as seen in the schematic trajectories in fig. 1. Subsequently, the evolution both in momentum and position has to be reversed. The second and third magnets first decelerate the particle and then impart momentum equal in magnitude to the momentum gained in the first stage, and finally the last magnet decelerates the particle, whereupon, if the conditions were ideal, we should find the partial beams reunited, and the particle in a pure spin state.
In our scheme, we exploit the correspondence between the Maxwell equations for the electric field and the static magnetic field, and introduce a fictitious magnetic "charge" potential φ, which obeys the Laplace equation and results in a physical magnetic field in full analogy to a static electric field. In this manner we obtain a configurable model of a magnetic field free from approximations. This is in contrast to other approaches to the Humpty-Dumpty problem, which tend to turn the problem of reversing the evolution into a technical difficulty depending solely upon the experimenter's skill, thereby possibly concealing difficulties of a more fundamental nature, while our aim is precisely to uncover those difficulties that arise once all simplifications are put aside. The scalar potential for the pair of magnetic line charges (see fig. 2) is where are the respective distances from the charged lines, and K 0 ( ) is the zeroth-order modified Bessel function of the second kind. On the Y axis, we have with b 0 (y) and λ(y) related to each other by where K 1 ( ) is the first-order modified Bessel function of the second kind. Upon expressing φ( r) in terms of b 0 (y), one confirms that φ( r) satisfies the Laplace equation (10) for s ± > 0 by an exercise in differentiation, for which the relations [39] K 0 (z) = −K 1 (z), are useful. After including the bias field B 0 e z and a strength parameter f , the magnetic field is Rather than specifying the line-charge density λ(y), we define the model by the choice for b 0 (y), essentially B = B z e z along the level trajectory of the atom that traverses the SGI. This makes it easier to ensure that the magnetic field decreases rapidly outside the SGA. Indeed, as we show below, one can model a magnetic field suitable for coherent beam recombination by an appropriate choice for b 0 (y), the strength parameter f , and the separation d = 2a between the two line charges.

Beam Recombination in the Calibrated Magnetic Field
The conventional treatment of an SGA assumes that the atom's deflection is caused by a linear B z ∝ z field, thus making the force F z ∝ ∂B z ∂z position independent, which greatly aids in subsequent beam recombination. However, a field like this is unphysical, it does not satisfy the Maxwell equations, and thus conclusions regarding the reversibility of the SGA's dynamics based on reasoning which neglects this issue are unconvincing as they risk missing quantum effects that might have an impact on the spin coherence. A physical field will necessarily have a non-linear dependence on the z coordinate and thus requires fine tuning to make the partial beams come back to zero separation.
In our scheme one can calibrate the field to result in a desired peak separation ∆Z by adjusting the field strength parameter f and the distance 2a between the top and bottom line charges. For b 0 (y) we choose a magnetic field profile composed of two gaussians with an adjustable width parameterized by α (as shown schematically in fig. 2 The parameter α determines the longitudinal stretch of the apparatus over which the particle is deflected, influencing the overall deflection, as well as the range of the fringing fields -the wider the gaussian peaks, the slower the decay of the magnetic field off the Y axis. As such, a judicious choice of α is needed, firstly to avoid a gaussian that is too narrow and expensive to resolve numerically, and secondly to avoid a gaussian too wide, which then makes it necessary to   fig. 1 in a simple way, without having to introduce a model of the magnetic field for each SG magnet separately. Note that this b 0 (y) is even in y and, therefore, we can replace the exponential factors in eq. (16) by their real parts.
The field felt by the particle along its trajectory is approximately B z ≈ b 0 (y), as long as the deflections are much smaller than the length scale of the field changes. Then, the vanishing of the integral makes the Larmor precession around the Z axis approximately equal to 0, with the spin pointing in the X direction at the output. Beam recombination along the Y axis is facilitated by energy conservation -changes in potential energy coming chiefly from the changes in B z ≈ b 0 (y) along the trajectory in the left half of the interferometer ( fig. 1) are reversed by the opposite changes in the right half (since b 0 (−y) = b 0 (y)). However, the choice b 0 (y) dictated by this reasoning is not a sufficient guarantee of coherent recombination. The real trajectory is not level and becomes minutely sensitive to the exact field, non-linear in its z coordinate dependence, in the apparatus. With the two remaining parameters at hand, the line charge separation 2a and the field strength f , we perform a bisection search for these parameters, constrained by the desired value of the peak separation ∆Z(T/2) ≡ ∆Z and with the objective of vanishing final spatial separation of the beams, ∆Z(T) ≤ 10 −6 δ.
This is done in a semiclassical manner: the Z axis quantum and classical dynamics are in good agreement, the classical simulation being much cheaper to perform than the quantum counterpart, thus making this optimization scheme feasible. We obtain a range of parameters a, f that result in spin-coherent recombination for different values of the for peak separation, ∆Z = 1δ, 2δ, . . . , 40δ.
With this set of optimized parameters {a, f } at hand (their values can be found in the appendix), each choice resulting in a suitable magnetic field and trajectories with peak separation ∆Z a, f , we can now proceed to simulating the quantum dynamics of the SGI.

Quantum Dynamical Simulation of the SGI
Using the optimized magnetic fields and the methods discussed in the preceding sections, we can now perform the wave-propagation simulation of the SGI. In the first step, we performed a quantum simulation of the beam's evolution in the optimized fields neglecting the microscopic quantum motions in the X, Y axes. The effective Hamiltonian for this case is a simplified onedimensional version of the full Hamiltonian of eq. (7): The simulated beam trajectories are presented in fig. 3. Regardless of the magnitude of the peak separation ∆Z = 1δ, 2δ, . . . , 40δ the beams are recombined coherently. Note that this happens only because the magnetic field, which depends on the z coordinate in a non-linear fashion (making the force F z position dependent), is finely calibrated for each ∆Z separately. The standard recipe of taking the "reverse" fields to reverse the evolution of the beam that underwent the splitting by SG magnets, therefore seems to work, as long as the fields are calibrated in an "exactly right" way. But is it really true? The standard argument assumes that the neglected X and Y movements are miniscule compared with the large deflection in the Z axis and, therefore, should not influence the overall analysis. However, since coherent beam recombination is concerned precisely with the relative microscopic realignment of the partial spin-up, spin-down beams, it is necessary to simulate the full dynamics.
We calculated the quantum dynamics of the original Hamiltonian of eq. (7), again varying the peak separation ∆Z and measuring the final spin coherence C(T). The peak separationspin coherence curve, together with the final relative displacements of the partial beams in y, z coordinates are presented in fig. 4. Although for very small separations up to ∆Z = 10δ the partial beams are recombined coherently-note that this does not mean that the evolution was fully reversed: the spreading of the wave packet, though inconsequential, has not been undoneonce the peak separation exceeds this value, the evolution rapidly becomes irreversible, with C(T) = 0.0 indicating complete spin decoherence by the time separation reaches ∆Z = 35δ. The cause of this loss of spin coherence is the growing microscopic separation in the Z and Y axes between the up and down spin components as shown in fig. 4b, which, although barely exceeding δ in Z and negligible in Y, is enough to completely prevent coherent beam recombination.

Discussion and Conclusion
An SGA entangles the spin degree of freedom of a spin-1/2 atom with its orbital degrees of freedom, so that the outcome of a position measurement informs us about the spin of the atom ("up" or "down" in the traverse direction that is probed). The transition from the initial not-entangled state to the final entangled state is unitary; it is not an irreversible event 4 that can be amplified and recorded. Since unitary processes compose a group, they are invertible in a mathematical sense. Is the entangling action of a SGA also invertible as a physical process -by another unitary transition, that is? 5 "Yes," answered Bohm [18] and later Wigner [19], and others parroted their wisdom. While Bohm acknowledged that the reversal would require a "fantastic" precision in controlling the apparatus and Wigner observed that such an experiment "would be difficult to perform," both took the reversibility for granted and did not elaborate on the matter. It is, however, unclear whether elementary unitary quantum processes are reversible -in marked contrast to macroscopic processes, which are known to be irreversible, as witnessed by the folk wisdom of the medieval Humpty-Dumpty rhyme, the laws of thermodynamics, the lessons of deterministic chaos, the ubiquitous ageing of living organisms, and other familiar phenomena.
In the particular situation of an SGA, the atom's spin and orbital degrees of freedom are coupled by the force resulting from the magnetic field gradient, and only this coupling is available for reversing the action of the SGA in a unitary fashion. We cannot undo the action of the unitary evolution operator exp −itH/h , with H as in eq. (3), by applying its inverse exp −it(−H)/h because −H is not a physical Hamiltonian -so much for mathematical reversibility. 6 Rather, we must carefully tailor the magnetic field to implement the disentangling reversal as well as we can.
We recall that Bohm's "fantastic" accuracy amounts to controlling the macroscopic devices with submicroscopic precision [20], and if such control is granted, Maxwell's equations prevent us from undoing the entanglement perfectly [21], even if we neglect the effects of fringing fields (as is the case in [21]). These theoretical works are here supplemented by a numerical study that fully accounts for the quantum dynamics with a magnetic field that obeys Maxwell's equations throughout the volume probed by the atom.
Thereby we assume that there are no uncontrolled imperfections. This is, of course, an over-idealization that lacks justification (as it ignores the lesson of [20]). But even with this stretch we find that the action of the SGA cannot be undone if the SGA serves its purpose, namely separates the atoms into well distinguishable partial beams -∆Z > 20δ, say.
In conclusion, the accurate wave-propagation results presented here quantitatively confirm that the microscopic quantum dynamical effects in a transversal SGI apparatus are enough to quickly destroy spin coherence beyond recovery. The deceptively simple nature of the Stern-Gerlach beam-splitting apparatus and the seemingly obvious proposal to reverse its evolution by employing "reverse" fields might lead one to assume the fundamental reversibility of its quantum evolution, as long as the apparatus is controlled sufficiently well. Such assumptions are, however, unwarranted.

Acknowledgments
We thank Jun Hao Hue for valuable discussions. The efforts by Tzyh Haur Yang and Ruiqi Ding, who conducted closely related undergraduate projects, are gratefully acknowledged.

Appendix: Particle properties and initial wave function
As in the original Stern-Gerlach experiment, we used a beam of silver atoms to simulate the splitting and recombination. The values of the various variables are given in table A1. In the lapse of time T, the width of the wave function would increase by a factor of if there were force-free motion, and this amount of spreading is of no consequence.

Appendix: Magnetic Field
The parameter α of the magnetic profile function b 0 (y) of eq. (19), depicted in fig. 2, was chosen to be The example parameters of the optimized fields, rounded to five decimal places, are displayed in table A2. In fig. A1 the field lines of the B z component of the field calibrated for ∆Z = 40δ (the largest simulated peak separation) are shown. The overall magnitude of the B y component of the generated field is actually larger than that the magnitude of B z . By adjusting B 0 we make sure that the resulting B z field component is the largest magnetic field component in the apparatus. The evaluation of the magnetic field integrals, the gradient of φ of eq. (16), across the wave function domain (the grid in the simulation for highest separation ∆Z = 40δ consists of 800 million points) was too time-consuming. Since the whole simulation was performed on a distributed computing cluster; we divided the domain into subdomains, with each subdomain being simulated on a separate processor. In each subdomain, then, the magnetic field B was approximated by a second-order Taylor approximation, centered at the weighted position of the wave function chunk in a particular domain, and there was no loss of accuracy from this.   Table A2. Optimized magnetic field parameters: beam separation ∆Z (cf. fig. 3), field-strength parameter f of eq. (18), and line-charge distance a (cf. fig. 2 The overall computational resources spent in the simulation of ∆Z = 40δ were as follows. The simulation was run on 40 cores in parallel, with a special parallelized FFT routine as the only step requiring inter-core communication. The computation time spent per core was 15 hours, a total of 600 hours. The breakdown of the average time spent in a single step of the time evolution loop: 30% spent on performing FFTs, 39% on potential propagators (costly because of field evaluations discussed in the preceding paragraph), 28% on observables calculations, and 3% on momentum propagators. The overall memory allocation per core was 3GB, a total of 120GB.  Figure A1. Lines of constant B z of the magnetic field calibrated for ∆Z = 40δ shown in the YZ plane. The sign of the field in a particular region is indicated by a + (B z > 0) and − (B z < 0). The nested contours increase in field strength as the field gets closer to the line charges located at a distance ±a at the top and bottom, according to table A2. Note that the width of the gaussian wave packet δ is small compared with the overall changes of the field. In the middle of the apparatus (Z = 0 line) the field vanishes exactly because the line charges of opposite sign are placed symmetrically across the XY plane.