Recovery Time of Matter Airy Beams using the Path Integral Quantum Trajectory Model

The Path Integral Quantum Trajectory (PIQTr) model is introduced as a new computational tool for the study of non-relativistic, quantum mechanical wave packets. We introduce the numerical algorithm and show that one of the primary advantages of the PIQTr model is its ability to efficiently handle heavy mass particles, while still resulting in a numerically exact answer. We then apply the model to the study of matter Airy beam particles and show that the recovery time of a damaged wave packet increases approximately with mass, but is independent of momentum, velocity, and wave packet width.


Introduction
There has been significant recent interest in the study and application of non-diffracting wave packets and twisted light and matter waves. With applications in fields such as micromanipulation, microscopy, and optical trapping, a fundamental understanding of the dynamics of these wave packets is crucial [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. One such type of non-diffracting wave packet is the Airy wave [16]. Unlike the more traditional Gaussian wave packets, Airy waves have many unique properties that offer features not seen with their Gaussian counterparts.
The first studies and applications of Airy wave packets have been in optical contexts [17][18][19][20][21][22][23][24][25][26][27][28], where they have been shown to be non-dispersive and exhibit force-free selfacceleration and self-healing properties. At first glance, some of these properties might seem too good to be true, however careful analysis shows that these features are real. For example, Airy beams are non-dispersive in the same sense that plane waves are non-dispersive. In their exact form, they are not square normalizable and extend throughout all space. They carry an infinite amount of energy [24] and its intensity remains constant throughout time. In reality, any Airy beam produced in the laboratory will be localized and therefore is square normalizable. This localization does cause the wave packet to diffract, although individual lobes within the beam do so at a slower rate than their Gaussian counterparts. The force-free acceleration seems to defy Ehrenfest's theorem, but while it is true that the individual spatial peak structures within the wave packet follow parabolic trajectories, the expectation value of a truncated Airy beam does not exhibit acceleration.
In addition to their propagation through free space, the propagation of optical Airy beams through nonlinear media has also been examined. In these media, Airy beams exhibit qualitatively different behavior, including features such as self-focusing [29][30][31][32]. Their stability over time is altered compared to similar wave packets in free space [33], and two interacting Airy beams have been shown produce spatial solitons [32].
While the vast majority of studies involving Airy wave packets have been in the optical context, recent experimental developments have demonstrated that electron Airy beams can now be produced in the lab [34]. These newly created matter Airy waves have many of the same properties as their optical counterparts and have the potential for many new and interesting applications, such as control and rotation of nanoparticles [6][7][8][9][10][11]19,35], improved resolution in electron microscopy [9,10,21,22], and the steering of electronic wave packets [34]. However, due to their recent discovery, matter Airy waves are much less understood.
Here we focus on the self-healing feature of Airy beams and define a recovery time for a damaged beam to recover its structure. The self-healing of an Airy beam occurs when a portion of the beam is perturbed or encounters an obstacle. Continued propagation of the beam transfers energy from the tail to the head of the beam, resulting in the main peak being restored, although at a decreased magnitude. We examine the propagation of damaged Airy waves in both free space and a nonlinear Kerr-type medium. We show that in both cases the recovery time increases as particle mass increases, and for wave packets in free space, it is independent of kinematical parameters. We also show that damaged Airy waves recover faster in a nonlinear medium than in free space.
With the recent creation of electron Airy beams, and the possibility for more massive Airy beams around the corner, it is essential to understand their unique properties at a fundamental level if the numerous proposed applications are to be achieved. To our knowledge, this work represents one of the first analyses of the self-healing properties of matter Airy beams in free space and in a nonlinear medium, and we hope that it spurs further theoretical and experimental developments in this area.
In order to study the dynamics of matter Airy beams, we introduce the Path Integral Quantum Trajectory (PIQTr) model as a new computational technique for finding the time evolution of particle wave functions and demonstrate the model's effectiveness in treating heavy particles. Because the PIQTr model does not rely on any classical, semi-classical, saddle point, or other approximation, it provides a highly accurate computational tool for the calculation of time-dependent wave functions and is especially useful for heavy particle dynamics.
Traditional approaches for studying the quantum mechanical motion and interactions of non-relativistic charged particles often focus on solving the Schrödinger equation. In the timedependent case, a variety of numerical methods are available, including finite difference methods, spectral representations, and finite element models. However, one of the challenges with both time-independent and time-dependent quantum mechanical methods is the computational difficulties encountered as the mass of the particle increases. In particular, methods based on basis expansions often require a large number of basis states, making the calculation difficult, if not impossible. Alternatively, matrix techniques can require the manipulation of large matrices, which is computationally challenging. In order to overcome these challenges, many methods employ classical or semi-classical treatments to make the calculation manageable.
One technique that is often overlooked when it comes to solving non-relativistic quantum mechanical dynamics is the Feynman path integral method, which has been widely used in high energy physics, quantum field theory, and statistical mechanics. However, it is not always a researcher's first instinct for fundamental physics problems [36][37][38][39][40][41][42][43][44]. Based on a Lagrangian approach, the formal path integral technique yields results completely equivalent to the Hamiltonian-based Schrödinger equation approach [45]. However, one of the practical difficulties associated with the method is the infinite sum over paths (or amplitudes associated with paths) required in the calculation. Our PIQTr model overcomes the need for this infinite sum by taking advantage of a heavy particle's small deBroglie wavelength in order to greatly reduce the computational requirements. The theoretical and numerical details of the PIQTr model are included in section 2, with a detailed analysis of the computational requirements presented in the appendices. Section 3 is dedicated to the discussion of matter Airy beam recovery time and in Section 4 a general summary and conclusion are given. Atomic units are used throughout unless otherwise stated.

PIQTr Model
The PIQTr model introduced here can be used to calculate time-dependent wave functions for a particle moving in any one-dimensional potential without approximation. This potential may or may not be time dependent. Fundamental quantum mechanics dictates that the time evolved wave function can be found by applying the propagator (also called the kernel or time evolution operator) to the initial state wave function where ψ( , ) is the initial state wave function at time ta, is the propagator and ψ( , ) is the new time evolved wave function at time tb.
The difficulty in the practical application of Eq. (1) is knowing the propagator. In the path integral method, the propagator is written as a sum over amplitudes for all trajectories for the particle to travel from its initial position to a final position. Each trajectory of the particle contributes an amplitude equal in magnitude, but different in phase is the classical action for the trajectory x(t) from xa to xb. It is given by where L is the Lagrangian of the system for a particle of mass m moving in a potential V(x,t). The constant A is the normalization constant. Combining Eqs. (1) and (2), the propagator can be expressed as where ( ) represents an integral over all trajectories. Inserting Eq. (5) into Eq. (1) yields At this point, the details of the PIQTr model become clearer by using a discrete formalism, with the added advantage of a direct connection to the computational method.
We begin by discretizing space and time, assuming for simplicity in this derivation uniform grids for both. Then, where x(t) represents an infinitesimal difference in paths.
Using Eq. (2), this becomes Consider now the evolution of the wave function through only one time step, such that ψ( , ) = ψ( , + Δ ). Then, assuming that the time step is infinitesimally small, there is only one path between two spatial points, which is a straight line. In this case, the propagator becomes 1 ( , , , + Δ ) = with S1step given by where L1step is the Lagrangian at time ta and position xa. In one dimension, the Lagrangian for a single time step is given by The normalization constant in Eq. (9) is given by 1 = √ 2 ℏΔ [1].
Combining Eqs. (8) - (11) yields the evolution of the initial state wave function by one infinitesimal time step In integral form, this is Equation (12) or (13) now allows for the iteration of the wave function by replacing Ψ( , ) with Ψ( , + Δ ). Successive applications will result in the time evolved wave function at any desired final time.
The derivation above is exact assuming that the time step Δt is small enough to approximate the limit that Δt -> 0. This requirement comes from the use of Eq. (10) and the assumption that for infinitesimal Δt there is only one straight-line path between two neighboring spatial points. A detailed discussion of the practical requirements for Δt can be found in Appendix A.
Brute force numerical evaluation of the integral in Eq. (13) is possible, but cumbersome due to the presence of K1step, which for fixed xb and Δt is an oscillatory function of xa. These oscillations increase as the magnitude of S1step increases, which occurs for small Δt, large mass, or large differences between xa and xb. However, it is possible to take advantage of the oscillatory nature of K1step when performing the integral. The principle of least action dictates that the classical path of the particle is the one that minimizes the action. Therefore, paths that deviate significantly from the classical path will have large values of the action, resulting in the single step propagator being a highly oscillatory function of xa. Assuming that the wave function is smooth over the wavelength of these rapid oscillations in the propagator, the integrand of Eq.
(13) is as likely to be positive as it is negative, and contributions to the sum of paths far from the classical path will effectively be zero [45]. This implies that the integration range for Eq. (13) does not need to encompass all of space, but rather a small region in which the single step propagator has minimum oscillations. The required integration range can be shown analytically to scale as δ ~ √ Δ , and numerical results agree with this to within a factor of 2 (see Appendix B for details).
The theoretical development above assumes constant temporal and spatial grids, and in practice, we use a constant time grid where Nt is the number of time grid points. However, the two spatial grids for the integration and the iterated wave function are not uniform and result from the Gaussian quadrature points used for the integration.
The runtime of the code scales linearly with Nt and as the product of the number of spatial integration points Nxint and the number of wave function grid points Nxwf. In general, the calculation times are quite manageable on a single processor, however parallelization of the integration is straightforward and serves to further improve runtime, scaling linearly with number of processors. A detailed analysis of the computational requirements of the PIQTr model can be found in Appendix C.

Results
As noted above, Airy beams were first found to be solutions to the free particle Schrödinger equation in 1979 [16], but only in the last decade have they received significant attention. Much of the focus has been within optical contexts [17][18][19][20][21][22][23][24][25][26][27][28] following the experimental and theoretical investigations of [17,28]. More recently, the production of electron Airy beams was demonstrated by passing electrons through a holographic mask that imparts a cubic phase [34]. These matter Airy beams exhibited the same unique properties of optical Airy beams, including curved trajectories, self-healing, and limited diffraction.

Free Particle Propagation
We begin with an analysis of free particle Airy wave packets. The discovery of Airy beams as non-diffracting particles was due to the relationship between the time-dependent free particle Schrödinger equation and the paraxial wave equation whose solution was known to be Airy functions [16]. The exact solution to Eq. (14) can be written as [16] ( , ) = [ where B is an arbitrary parameter (taken to be 1 here). From the t 2 term in the Airy function, it is easy to see that this wave will exhibit an acceleration of 1 2 2 [16]. In our investigations, we use an initial state, exponentially truncated Airy beam with the leading peak centered at x = 0 and general form [18] The truncation parameter determines the width of the initial Airy wave packet. A larger value of results in a more narrow initial wave packet. Because the acceleration of the Airy beam peaks is known, the classically predicted trajectory is given by Clearly, as the mass of the particle increases, the acceleration is reduced. Figure 1 shows the time evolution of free particle Airy beams with truncation parameter = 0.1, initial velocity v0 = 1 and masses m = 1, 10, and 100 a.u. With time on the vertical axis and space on the horizontal axis, the color scale represents the probability density | ( , )| 2 .
The decreasing acceleration with increasing mass is apparent with the trajectories for heavier particles having less curvature. It can also be seen that the Airy wave packet maintains its shape and magnitude for a finite time. Over time, the Airy wave packet spreads, leading to a decrease in the truncation parameter for > 0. Calculations show that the truncation parameter is halved at a time that is directly proportional to mass. Therefore, the stability of an Airy beam's shape over time is improved for more massive particles. The primary feature of the matter Airy beams in which we are interested is their selfhealing nature. When the Airy beam is perturbed or damaged due to an obstacle, the beam will recover its initial structure after some time. This has been demonstrated extensively with optical Airy beams [17,18,23,24], and shown to hold true with electron Airy beams [34]. Given the recent intense interest in optical Airy beams and the creation of electron Airy beams, it is worthwhile to determine how the self-healing properties change with particle mass. In particular, because any practical application of matter Airy beams will involve their interaction with matter, it is important to quantify the time it takes for the wave packet to recover its shape following some distortion or damage to the original wave packet. For example, if an Airy wave packet is perturbed through a scattering event or truncated by a collimating slit, it will be useful to know if the wave packet will recover its shape before a second interaction might occur. In particular, we study Airy beams whose initial leading peak (in the direction of travel) has been blocked by an obstacle, as might be encountered with a collimating slit for a scattering experiment. For comparison, Fig. 2 shows initial state damaged and undamaged Airy beams.
Both are normalized such that ∫ | ( , 0)| 2 = 1. In order to quantify the time it takes an Airy beam to self-heal, we introduce a definition of 'recovery time', which we define as the time it takes for the location of the leading peak maximum in the damaged Airy beam to reach the predicted location of the leading peak maximum in an undamaged Airy beam. In other words, the time for the trajectory of the damaged Airy beam's leading peak to reach the predicted path of an undamaged Airy beam's leading peak is defined as the recovery time.
This definition of recovery time presents a practical challenge because the damaged Airy beam's leading peak approaches the predicted location gradually in an asymptotic-like fashion.
Thus, a specific time for recovery is difficult to define. However, a general trend of how recovery time scales with mass can be found. We determine the recovery time by taking the ratio of the predicted leading peak location to the leading damaged peak's location. When this ratio is within 1% of unity, we consider the damaged Airy beam to be fully recovered. Fig. 3 shows a sample plot of the ratio of leading peak locations, as well as the leading peak trajectory of the damaged Airy beam and the corresponding undamaged trajectory (m = 1, p0 = 1, = 0.1).
The discontinuity in the ratio occurs at the time where the damaged Airy beam's leading peak location changes from a negative to a positive value. The undamaged leading peak's position is always positive. At long times, the peak location ratio approaches unity as the damaged Airy beam returns to its original trajectory and remains one for further propagation. This indicates that the Airy beam has recovered and will follow its original trajectory indefinitely. The blue arrows shown in Fig. 3 point to the recovery time when the peak location ratio is within 1% of unity. . . ). Figure 4 shows the free particle time evolution of three damaged Airy beams with truncation parameter = 0.1. Each of the Airy beams has been damaged by removing the leading peak (see Fig 2). Results are shown for particles with (m,p0) = (1,1), (10,10), and (1,10).

Figs 4 (a) and (b)
show results for particles with the same velocity, but different masses. A comparison of these two cases shows that the recovery time varies with mass, despite the same particle velocity . Figs 4 (a) and (c) show results for particles with the same mass, but different velocities and momenta, and it can be seen that the recovery time does not change as the initial velocity changes. In addition to varying the mass, momentum, and velocity, we also performed calculations for Airy beams with a larger truncation parameter and found that the recovery time was approximately independent of width, as shown in Table 1. As can be seen from Figs 3 and 4, the process of self-healing evolves such that the leading peak in the damaged wave function gradually returns to the trajectory of the leading peak for the undamaged wave function (shown as the pink line with circles). Quantitative results for recovery times are shown in Table 1, and it can be seen that for free particle propagation the recovery time scales approximately linearly as four times the mass of the particle. In general, all of our testing indicates that the recovery time is primarily dependent upon the mass of the particle and is independent of kinematical parameters such as momentum, velocity, and truncation. The approximately linear scaling of the recovery time with mass and the insensitivity to momentum, velocity, and initial width implies that more massive Airy beams, such as those created with muons or protons, will take a very long time to recover. This could have important implications on any applications that might rely on the self-healing properties of Airy beams.  Table 1 Recovery time of damaged Airy beams with initial velocity of 1 a.u., but different masses and truncation parameters.

Nonlinear Propagation
The studies of optical Airy beams in free space have been expanded to their propagation and interaction in nonlinear media, such as a Kerr medium in which the presence of an electric field alters the medium's index of refraction [32,33,[46][47][48]. Propagation of a single optical Airy beam in a focusing Kerr medium results in an enhancement of the second lobe intensity [33], but the leading lobe follows the same trajectory as in free space. The nonlinearity can also affect the stability of the Airy beam [48][49][50] and interacting beams have been shown to generate spatial solitons [32]. Additional features such as self-focusing or trapping of optical Airy beams in nonlinear media [51] have also been observed.
In general, matter Airy beams have been studied much less than those in optical contexts, and the same is true for their propagation in nonlinear media. While the Kerr effect is an optical response, there is evidence for the propagation of matter waves in nonlinear media. For example, the nonlinear Schrödinger equation with a cubic interaction term (as in the Kerr effect) can be used to describe some Bose Einstein condensates [52,53].
Here we investigate matter Airy beams moving in a nonlinear focusing Kerr-type medium where the wave function evolves according to the nonlinear Schrödinger equation where is given by | ( )| 2 and g is a nonlinearity constant. In an optics context, represents a variable index of refraction, while for matter Airy beams, is a mean field term that is proportional to particle density [53].
As in the optics case, the propagation of matter Airy beams in the Kerr-type medium is qualitatively different than that of free particles. Figure 5 shows a comparison of undamaged free space Airy waves to those in the Kerry-type medium for m =1 and m = 10. The matter Airy wave in the nonlinear medium shows an enhancement in the secondary peaks at small times for both masses. As the beams are allowed to evolve, the leading peak in the nonlinear medium spreads and loses intensity over time, unlike the wave packet in free space. Also, in the Kerrtype medium, the leading peak maximum follows a similar trajectory to that of an Airy wave packet in free space, but lags slightly behind (see Figs 5 and 6). This decreased acceleration could be due to the self-focusing that occurs in nonlinear media [48], where energy intensifies each lobe, in addition to accelerating it. Also, some interference structures can be observed in the tail of the wave function at longer times. These structures result from negative momentum components that enter into the wave function as a result of the nonlinear medium. Next, we turn to the calculation of recovery time for the Airy beam in a Kerr-like medium. We use the same definition of recovery time that was introduced above and note that the dynamics of the evolving damaged Airy beam in the nonlinear medium are very similar to those of the damaged Airy beam in free space. However, the evolution of the undamaged Airy beam in the nonlinear medium is different from its free space counterpart. This leads to differences in the recovery times. In particular, the damaged Airy beam no longer asymptotically approaches the undamaged beam trajectory in the Kerr-type medium. Instead, the damaged Airy beam leading peak trajectory approaches and then crosses the undamaged trajectory, asymptotically approaching the undamaged free space trajectory, as seen in Fig 6. The results in Fig 6 are shown for m = 1 and p0 = 1, but results for other mass particles show similar behavior. Because the undamaged beam accelerates more slowly in the nonlinear medium, the damaged Airy beam recovers more quickly than it does in free space, with recovery times 25 to 50% faster (see Table 2).  Table 2 Recovery time of damaged Airy beams with initial velocity of 1 a.u., = 0.1, and different masses for propagation in the nonlinear Kerr-type medium and free space.

Conclusion
In summary, we examined the self-healing property of damaged Airy beams in free space and a nonlinear Kerr-type medium. Damage was induced by removing the leading peak to simulate an encounter with an obstacle. A recovery time for the wave packet was defined and shown to increase approximately linearly with increasing particle mass in free space.
Additionally, in free space, this recovery time was independent of both velocity, momentum, and width of the wave packet. In the nonlinear medium, the recovery time also increased with particle mass, but not linearly. A comparison between recovery time in free space and recovery time in the nonlinear medium showed that the Airy beam recovers more quickly in the nonlinear medium. This decreased recovery time could be due to the self-focusing aspects of Airy beam dynamics in nonlinear media. The increase in recovery time with mass will have important implications in any future experiments involving massive Airy beams and indicates that even in nonlinear media, the heavy Airy beams will be unlikely to self-heal during the time frame of an experiment.
The time-dependent Airy beam wave functions were calculated with our newly introduced Path Integral Quantum Trajectory (PIQTr) model and the numerical algorithm of the method was introduced. By taking advantage of the smaller deBroglie wavelength of heavy particles, the PIQTr model scales favorably with increasing mass. This feature will be useful for future applications to heavy particle dynamics and lays the foundation for the PIQTr model being one of the first that can be applied to both heavy and light particles with similar computational requirements.

Appendix A
The theoretical derivation of the PIQTr model requires that the time step for iterations of the wave function be infinitesimally small such that there is only one linear path between two spatial points. Then, the propagator can be written as Eq. (9). The corresponding numerical requirement is that Δ is sufficiently small that the propagated wave function will not change as Δt decreases. Experience shows that if Δt is too large, the wave function will propagate too quickly. However, if it is too small, the computation becomes intensive. The key to achieving accurate results is to keep Δt as small as possible without resulting in an unnecessary computational burden. Unfortunately, there is no way to analytically approximate Δt because it will depend on the potential chosen, as well as the mass.
We can gain, however, insight into the requirements of Δ by studying the familiar Gaussian wave packet moving under the influence of a constant force. In this case, the analytical result is well-known, which allows for detailed testing of the effect of Δt on the time propagation of the wave function. Other systems will have similar, but not identical, requirements. The force is chosen to be 2 a.u. and the mass is chosen to be 1 a.u. The initial state wave function is a Gaussian with σ = 1, no initial momentum (p0 = 0), and initial position x0 = 0 a.u.     From Fig. A2, it is easy to see that as the time step decreases, the agreement between the exact answer and the numerical result improves. When Δt is too large, the K1step peaks are shifted to larger xa from Kexact. Also, as Δt increases, the peaks of Kexact shift away from zero.
Therefore, a Δt value that is too large effectively gives peaks in locations similar to those of Kexact at a later time, which results in an increase in velocity. Finally, in Fig. A2, one also sees that as the time step increases, fewer oscillations are observed in the propagator over the same spatial distance. This is easily predicted from Eqs. (10) and (11), where the action is seen to vary inversely with Δt in the first term and linearly with Δt in the second term. Thus, as Δt increases, there are fewer oscillations over the same distance, indicating that a particle may move more slowly in order to arrive at the same xb point at a later time.

Appendix B
In its exact form, the PIQTr model requires a numerical integration over all space, however in practice the finite deBroglie wavelength of a particle allows for this integration region to be reduced. To determine the required numerical integration region, note that K1step is centered around xa = xb. When xa is sufficiently far from xb, the integrand becomes highly oscillatory and the positive and negative contributions to the integral will cancel. Therefore, for |xa -xb| > δ, there is no need to perform the integration and Eq. This results in a significant computational advantage due to the much smaller integration range.
A semi-empirical estimate for δ can be found from examining how the wavelength of the single step propagator changes with increasing |xa -xb|. Consider the free particle single step propagator in which V(xa,ta) = 0 in Eq. (11). The de Broglie wavelength can be written as We assume that it is safe to neglect the contributions to the integral when the change in wavelength as a function of xa is less than some value β. Therefore, The (xbxa) 2 term in the denominator is simply δ 2 and therefore, From this expression, it is apparent that δ ~ √ Δ for fixed β. Recall that Δt is found numerically from the approximation of the action integral and β can be assumed to be a constant independent of mass that determines when K1step has approximately constant wavelength oscillations. From numerical testing of a free particle with mass 1, we found that δ = 8 and Δt = 0.05 are sufficient for a numerically exact answer. This results in β = 0.005. Therefore, as mass increases, δ decreases such that the range of integration also decreases. Also, as Δt decreases, δ decreases and a smaller integration range may be used. It is in this manner that the PIQTr method is able to take advantage of a decreasing de Broglie wavelength for increasing mass. Table B1 shows the predicted values of δ from Eq. (B4) compared with the values found from numerical testing.
While not exact, Eq. (B4) serves as a guide to estimate δ, generally to within a factor of 2.
Mass δ estimate from (17) δ found numerically 10 3 3 100 1.1 2 1000 0.79 1.5 10000 0.35 0.5 Table B1 Comparison of necessary integration range as predicted from Eq. (B4) with results obtained from numerical testing. Note that a mass of 1 was used to find β, which was in turn used to calculate δ, and therefore the two values must be identical. It was therefore excluded from the table.

Appendix C
One of the major challenges with the numerical calculation of heavy particle wave functions is the numerical and computational difficulty associated with increasing mass. It is therefore important that we understand how the computational requirements scale with increasing mass. We show in Table C1  From Table C1 it can be seen that the runtime per iteration changes by about a factor of 30 as the mass of the particle increases by 4 orders of magnitude. The primary parameters that influence the runtime are Nxwf and Nxint. However, the minimum value for these parameters is coupled to the mass, Δt, and the range of integration δ. A larger mass or smaller Δt leads to greater oscillations in K1step, which will require a greater density of Nxint and Nxwf points.
However, a smaller δ for the same number of Nxint and Nxwf points leads to a greater point density.
Therefore, as mass increases or Δt decreases, we are able to keep Nxint and Nxwf at manageable values due to the decreasing δ required. We note that while it is advantageous for the runtime per iteration to remain relatively stable, a particle with larger mass will require a larger final time to reach the same final position. Therefore, the total computational time for the larger mass particles is increased over the smaller mass particles. However, the required minimum value for Δt increases as mass increases, which helps keep runtimes down, even for large mass particles.