#
Possible potentials responsible for stable circular relativistic orbits

Prashant Kumar†, Kaushik Bhattacharya‡
Department of Physics, Indian Institute of Technology, Kanpur
Kanpur 208016, India email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Bertrand’s theorem in classical mechanics of the central force fields attracts
us because of its predictive power. It categorically proves that there can
only be two types of forces which can produce stable, circular orbits. In the
present article an attempt has been made to generalize Bertrand’s theorem to
the central force problem of relativistic systems. The stability criterion for
potentials which can produce stable, circular orbits in the relativistic
central force problem has been deduced and a general solution of it is
presented in the article. It is seen that the inverse square law passes the
relativistic test but the kind of force required for simple harmonic motion
does not. Special relativistic effects do not allow stable, circular orbits in
presence of a force which is proportional to the negative of the displacement
of the particle from the potential center.

##  1 Introduction

The central force problem in non-relativistic classical mechanics is one of
the most useful topics in physics. Closely linked with the central force
problem is the Keplerian orbit theory which is a cornerstone for understanding
planetary motions in the solar system or motion of electrons near the nucleus.
In classical mechanics there is an important theorem called the Bertrand’s
theorem which proposes that there can only be two types of central potentials,
the Coulomb type and the simple harmonic type, which can produce stable,
circular orbits for particles moving around the potential source. A good
presentation of the Bertrand’s theorem can be found in Ref. [1]. The present
article tries to generalize the results of Bertrand’s theorem when the
orbiting particle can have relativistic velocities.

In this article we first set up the relativistic orbit equation for a particle
in a central potential presumed to be dependent on the radial coordinate only.
The relativistic central force orbits were previously studied in Refs. [2, 3,
4, 5]. A brief description of the central force problem in a relativistic
setting in a Coulomb potential was presented in the book on classical theory
of fields by Landau and Lifshitz [6]. Before one starts the main analysis
about the stability of orbits of relativistic particles in a central force
potential it is better to specify the assumptions one makes in arriving at
definite results. In the present article we use the same assumptions and the
approximations as utilized by Boyer in Ref. [2] and Landau in Ref. [6]. In the
specific references cited above, none of them present a Lorentz covariant
treatment of the relativistic central force problem. The main reason being
that all of them assumes a central potential $V(r)$ where $r=|{\bf r}|$ is the
distance between the source and the orbiting particle. The form of the
potential only depends on the position coordinates of the orbiting particle.
The form of $V(r)$ is not Lorentz covariant. In such cases the results of the
whole analysis is valid in a particular frame where the origin of the
coordinate system coincides with the potential center.

The references cited above assumes the particle which produces the potential
$V(r)$ to be static in the specific coordinate system utilized by the
observer. If the source of the potential does not have any velocity then the
retarded nature of the interactions, owing to the finite velocity of light,
does not complicate the calculation of the orbit of the relativistic particle.
A specific example will make the point clear. In classical electrodynamics if
the source of the Coulomb potential $V(r)$ moves with a velocity ${\bf v}_{s}$
and the orbiting particle has a velocity ${\bf v}$ then the potential $V(r)$
gets a relativistic correction. The magnitude of the lowest order relativistic
correction to the Coulomb potential was calculated by Darwin in 1920 and it
looks like

$\frac{V(r)}{2c^{2}}\left[{\bf v}_{s}\cdot{\bf v}+\frac{({\bf v}_{s}\cdot{\bf
r})({\bf v}\cdot{\bf r})}{r^{2}}\right]\,.$

For a better understanding of the Darwin correction one can look at Ref. [7].
In our case ${\bf v}_{s}=0$ and consequently there will be no relativistic
modification of $V(r)$. More over we do not consider any general relativistic
effects due to $V(r)$ into account. We briefly comment on the general
relativistic generalization of the central force problem in section 4. In the
present article the background space-time is assumed to be flat.

In the article it will be shown that the stability condition of the perturbed
orbits around a stable circular orbit gives rise to a non-linear differential
equation for the central potential. The Newtonian or the Coulomb potential
satisfies the resulting differential equation with some restrictions on the
possible value of the angular momentum of the orbiting particle. Except the
Newtonian potential solution we present a more general solution of the
differential equation for the potential which can give rise to stable,
circular orbit for relativistic particles. This solution gives rise to a force
which is not common in physics except its Newtonian inverse square law limit.
The equation of the orbit of a relativistic particle in such a non-trivial
force shows that the orbit will precess and the precession angle can be
calculated.

Unlike the non-relativistic case, in the relativistic case there exist no
radial effective potential minimizing which we can obtain the radius of a
circular orbit. In the relativistic case a first order perturbation from a
circular orbit is enough to determine the stability criterion of the orbit. In
the non-relativistic case one uses higher order perturbations from a circular
orbit to specify the form of the potential. In the relativistic case the
general solution of the form of the potential from first order perturbation
from circular orbit is such that all higher order corrections becomes
irrelevant. As a consequence of this fact the general form of the potential
which can produce stable, circular orbits for relativistic particles contains
more parameters than the corresponding expressions of non-relativistic
potentials.

The material in the article is presented in the following manner. The second
section sets the conventions and derives the orbit equation of a relativistic
particle in a central orbit. Section 3 generalizes the Bertrand’s theorem for
the relativistic case. In this section the stability condition for the
circular orbits will be interpreted as a non-linear differential equation for
the potential. The solutions of the non-linear stability equation will also be
derived in section 3. In section 4 the connection of the present work with
some related works which were existing in the literature are discussed. This
section gives a wider view for the readers who really want to understand the
stability of orbits in special relativity and general relativity. The last
section 5 summarizes the important points presented in the article.

##  2 The orbit equation

In this section we derive the orbit equation of the relativistic particle in
presence of a potential $V(r)$ which is purely a function of the radial
coordinate. The Lagrangian of a relativistic particle of mass $m$ in presence
of a continuous radial potential $V(r)$ is

$\displaystyle{\mathcal{L}}$ $\displaystyle=$ $\displaystyle-
mc^{2}\sqrt{1-{v^{2}}/{c^{2}}}-V(r)\,,$ (1)

where the velocity of the particle ${\bf v}$ in plane polar coordinates is
given as

${\bf v}=\dot{r}\hat{e}_{r}+r\dot{\theta}\hat{e}_{\theta}\,,$

where $\hat{e}_{r}$ , $\hat{e}_{\theta}$ are the mutually orthogonal unit
vectors along the radial and the angular directions. The form of the
Lagrangian in Eq. (1) immediately shows that the angular momentum

$\displaystyle
L=\frac{\partial{\mathcal{L}}}{\partial\dot{\theta}}=mr^{2}\gamma\dot{\theta}\,,$
(2)

is a constant, where

$\gamma=\frac{1}{\sqrt{1-{v^{2}}/{c^{2}}}}\,.$

The total energy $E$ of the particle in presence of the potential $V(r)$ is

$\displaystyle E=mc^{2}\gamma+V(r)\,.$ (3)

Although from the definition of $\gamma$ it looks like that it is a function
of $r$, $\dot{r}$ and $\dot{\theta}$ but it can be shown that in a central
force field $\gamma$ is only a function of the radial coordinate $r$. The
reason for such behavior of $\gamma$ can be understood from the following
reason. As energy and angular momentum are constant functions of $r$,
$\dot{r}$ and $\dot{\theta}$ we can use the conservation conditions of $E$ and
$L$ to re-express $\dot{r}$ and $\dot{\theta}$ as functions of $r$, $E$ and
$L$. As $E$ and $L$ are constants so in a central force field $\dot{r}$ and
$\dot{\theta}$ are functions of $r$ alone. Consequently $\gamma$ is only a
function of $r$. In special relativity the energy of the particle in a central
potential can also be written as

$\displaystyle E=\sqrt{p^{2}c^{2}+m^{2}c^{4}}+V(r)\,,$ (4)

where $p=|{\bf p}|$,

$\displaystyle{\bf p}=m\gamma{\bf v}$ $\displaystyle=$ $\displaystyle
m\gamma(\dot{r}\hat{e}_{r}+r\dot{\theta}\hat{e}_{\theta})$ (5)
$\displaystyle=$ $\displaystyle
p_{r}\hat{e}_{r}+p_{\theta}\hat{e}_{\theta}\,,$

and

$p_{r}=m\gamma\dot{r}\,,\,\,\,\,\,\,p_{\theta}=m\gamma
r\dot{\theta}=\frac{L}{r}\,.$

As because $({p_{r}}/{p_{\theta}})=({\dot{r}}/{r\dot{\theta}})$, we have

$p_{r}=\frac{L}{r^{2}}\frac{dr}{d\theta}\,.$

With the above information on the various momentum components we can now
rewrite Eq. (4) as

$\displaystyle(E-V)^{2}=\left(\frac{L}{r^{2}}\frac{dr}{d\theta}\right)^{2}c^{2}+\frac{L^{2}c^{2}}{r^{2}}+m^{2}c^{4}\,.$
(6)

Instead of $r$ we use the variable $u=\frac{1}{r}$ in terms of which Eq. (6)
becomes

$(E-V)^{2}=L^{2}c^{2}\left(\frac{du}{d\theta}\right)^{2}+u^{2}L^{2}c^{2}+m^{2}c^{4}\,.$

If we differentiate the last equation with respect to $\theta$ and then divide
the resulting equation by $du/d\theta$ we obtain the desired equation of the
orbit of a particle of mass $m$ possessing momentum ${\bf p}$ moving in the
presence of a general central potential $V(r)$ as

$\displaystyle\frac{d^{2}u}{d\theta^{2}}+u=\frac{(V-E)}{L^{2}c^{2}}\frac{dV}{du}\,.$
(7)

Using Eq. (3) we can rewrite the above equation in the form

$\displaystyle\frac{d^{2}u}{d\theta^{2}}+u=-\frac{m\gamma}{L^{2}}\frac{dV}{du}\,.$
(8)

Writing $L=\gamma\ell$, where $\ell=mr^{2}\dot{\theta}$ is the non-
relativistic angular momentum, the above equation in the non-relativistic
limit ($\gamma\to 1$) transforms exactly to the form we get in a conventional
non-relativistic treatment of the problem as given in Ref. [1].

##  3 Circular, Stable closed orbits

Lets define

$\displaystyle J(u)\equiv\frac{(V-E)}{L^{2}c^{2}}\frac{dV}{du}\,.$ (9)

Suppose Eq. (7) admits a circular orbit of radius $r_{0}=1/u_{0}$. For small
perturbations around this circular orbit we can Taylor expand $J(u)$ around
$u_{0}$. Keeping up to first order terms in the perturbation of $u$ we get

$\displaystyle J(u)=J(u_{0})+(u-u_{0})\left(\frac{dJ}{du}\right)_{u_{0}}\,.$
(10)

Noting that $J(u_{0})=u_{0}$ for the circular orbit, we can now write Eq. (7)
as

$\displaystyle\frac{d^{2}u}{d\theta^{2}}+(u-u_{0})=(u-u_{0})\left(\frac{dJ}{du}\right)_{u_{0}}\,.$

If we define $x\equiv u-u_{0}$ then the above equation can be written as

$\displaystyle\frac{d^{2}x}{d\theta^{2}}+\zeta^{2}x=0\,,$ (11)

where $\zeta^{2}$ is defined as

$\displaystyle\zeta^{2}\equiv 1-\left(\frac{dJ}{du}\right)_{u_{0}}\,.$ (12)

From Eq. (11) it is clear that if the orbit of the relativistic particle in a
general central potential has to be stable then $\zeta^{2}>0$ and if the orbit
has to be closed then $\zeta$ must be a rational number.

###  3.1 A differential equation for the potential $V(r)$ producing stable and
closed circular orbits

The rational number $\zeta$ as predicted, in Eq. (12), from the stability
criterion of closed circular orbits in the central force problem is an
interesting input in the theory. The interesting property about this rational
number is that it is a constant and so it does not depend on the details of
the orbit which one tries to perturb. The reason for the constancy of $\zeta$
is the following. For any circular orbit with radius $r_{0}$ a specific
$\zeta$ specifies the number of undulations of the perturbed orbit. If $\zeta$
is a rational number then the number of undulations of the perturbed orbit
will be such that they form a closed geometrical structure. Now suppose one
takes another circular orbit of radius $r_{0}+\delta r$ where $\delta r\ll
r_{0}$. If $\zeta$ has a different value on this orbit then the number of
undulations due to a perturbation will be different. In the limit $\delta r\to
0$ in a continuous manner the two unperturbed circular orbits tends to each
other but the number of undulations on the circular orbits will not match as
$\zeta$ is not a continuous variable but can only have discrete rational
values. Consequently the number of cycles of the perturbations will change
discontinuously with radius and the perturbed orbits cannot be closed at this
discontinuity. As we are only interested in stable, closed orbits we can
conclude that $\zeta$ must be a constant and not change discretely with $r$.
The discussion on the constancy of $\zeta$ as given above closely follows the
analysis given in Ref. [1] where the author gives a nice discussion on the
role of $\zeta$ in the case of non-relativistic orbits.

As $\zeta$ is a constant and must not depend upon the choice of $u_{0}$ or $x$
one can interpret Eq. (12) as an independent differential equation by itself,

$\displaystyle 1-\left(\frac{dJ}{du}\right)=\zeta^{2}\,.$ (13)

whose solutions would give us information about the general form of the
central potential $V(r)$. Using Eq. (9) we can write the last equation as

$\displaystyle(V-E)\frac{d^{2}V}{du^{2}}+\left(\frac{dV}{du}\right)^{2}=L^{2}c^{2}(1-\zeta^{2})\,,$
(14)

which is a non-linear second order differential equation. The right hand side
of the above equation is a constant which can be written as

$\displaystyle d=L^{2}c^{2}(1-\zeta^{2})\,.$ (15)

Eq. (14) admits multiple solutions for $V$. The constant $E$ is the total
energy of the particle.

It is interesting to note that the differential equation for the potential
stemming from the stability of closed, circular orbits in the relativistic
case does not have a non-relativistic analogue. Although the orbit equation
Eq. (8) has a proper non-relativistic limit the same cannot be said about Eq.
(14). The reason for such behavior can be seen clearly if we rewrite Eq. (14)
in a slightly different way. From the expression of the energy of the particle
in the central force field $\gamma$ can always be written as $(E-V)/mc^{2}$.
As the total energy is a constant in the present case we must have
$d\gamma/du=-(1/mc^{2})dV/du$. Consequently Eq. (14) can also be written as

$\displaystyle\gamma\frac{d^{2}\gamma}{du^{2}}+\left(\frac{d\gamma}{du}\right)^{2}=\frac{L^{2}(1-\zeta^{2})}{m^{2}c^{2}}\,,$
(16)

which gives a differential equation of $\gamma$. The equation above obviously
does not have a well defined non-relativistic limit. The relativistic
stability condition produces an ill-defined non-relativistic limit due to the
fact that in the relativistic case $J(u)$ as given in Eq. (9) depends upon the
velocity of the orbiting particle111The $V-E$ in $J(u)$ is proportional to
$\gamma$ which depended upon the velocity of the particle.. In the non-
relativistic $J(u)$ was purely a function of the radial coordinate of the
orbiting particle. A perturbation from the circular orbit in the relativistic
case consists of two kinds of perturbations. One is related to the change in
position of the particle from its previous orbit and the other is the change
in velocity from the velocity it had previously on the circular orbit. In the
non-relativistic case only a radial perturbation from the circular orbit fixes
the shape of the stability condition. As because the stability condition of
the orbit depends upon velocity of the relativistic particle and the
corresponding non-relativistic stability condition does not depend upon the
velocity of the particle, the non-relativistic limit of Eq. (14) or Eq. (16)
is not well defined.

In the case of non-relativistic motion we know that the inverse square law
potential and the simple harmonic potential has the capability to produce
stable, closed circular orbits. In the present case to get the forms of the
potentials which can produce stable, closed orbits we have to solve Eq. (14).
As it is a non-trivial equation we will first try to see whether the the
potentials which produced stable, closed orbits in the non-relativistic regime
still satisfy Eq. (14). Let us try to see whether any power law solution of
the form

$\displaystyle V(u)=-\alpha u^{\tau}\,,$ (17)

where $\alpha>0$ satisfies Eq. (14). In the above equation $\alpha$ and $\tau$
are constants. If we substitute the above form of the potential in Eq. (14) we
get

$u^{2(\tau-1)}\left\\{\alpha^{2}\tau(\tau-1)+\alpha^{2}\tau^{2}\right\\}-u^{\tau-2}\left\\{E\alpha\tau(\tau-1)\right\\}=d\,.$

This directly shows that the above relation can be valid for any $u$ only if
$\tau=1$, when $\alpha^{2}=d$ or

$\displaystyle L=\frac{\alpha}{c\sqrt{1-\zeta^{2}}}\,,$ (18)

on using Eq. (15). In this case we see that choosing $\tau=1$ in Eq. (17) we
get the Coulomb or Newtonian potential. The last equation shows that for
stable, circular orbits the particles angular momentum must satisfy some
condition. Eq. (18) implies that $\zeta^{2}<1$, and as $\zeta^{2}>0$ for a
stable orbit, we have

$\displaystyle 0<\zeta^{2}<1\,.$ (19)

The above equation gives

$\displaystyle L>\frac{\alpha}{c}\,,$ (20)

giving a lower bound on the angular momentum of the orbiting particle. This
lower bound of the orbiting particle was previously obtained in a different
way by T. H. Boyer in Ref. [2]. It must be noted here that except $\tau=1$ no
other values of $\tau$ are allowed in the potential which can produce stable
circular orbits of relativistic particles . In non-relativistic mechanics we
do also have the harmonic-oscillator potential corresponding to $\tau=-2$ and
$\alpha<0$ in Eq. (17), but interestingly relativistic effects forbid this
value of $\tau$.

###  3.2 The general solution of the differential equation for the potential
and the nature of orbits

We can find out the general form of the force which can produce stable,
circular relativistic orbits. Noticing that the left hand side of Eq. (14) can
also be written as

$\frac{d^{2}}{du^{2}}\left[\frac{(V-E)^{2}}{2}\right]\,,$

it can be easily shown that

$\displaystyle V(r)-E=-\sqrt{d\left(b+\frac{1}{r}\right)^{2}+a}\,,$ (21)

satisfies Eq. (14) where $d$ is as given in Eq. (15) and $b$ and $a$ are two
other dimensional, integration constants. For an attractive force $b>0$ and
$d>0$ but $a$ can have any sign. If we assume that as $r\to\infty$, $V(r)\to
0$ then we get a relation between the constants $d$, $b$ and $a$ as

$\displaystyle E=\sqrt{db^{2}+a}\,.$ (22)

From Eq. (21) we get the force acting on the particle,

${\bf F}=-\nabla V(r)\,,$

as

$\displaystyle{\bf
F}=-\frac{d\left(b+\frac{1}{r}\right)}{r^{2}\sqrt{d\left(b+\frac{1}{r}\right)^{2}+a}}\,\,\hat{r}\,,$
(23)

From the form of the force and Eq. (22) we immediately see that if $a=0$ we
have $b=E/\sqrt{d}$ and we get back the Newtonian or the Coulombic potential.
From the form of the potential as written in Eq. (17) we can furthermore
identify $\alpha=\sqrt{d}$ and consequently when $a=0$ we have $b=E/\alpha$.
If $a\neq 0$ then the form of the force is non-trivial. The form of the force
as given in Eq. (23) cannot be reduced to the harmonic oscillator force in any
limits of the constants. This shows that special relativistic effects do not
allow stable circular orbits in the presence of a force which is proportional
to the negative of the displacement vector. Although the force expression in
Eq. (23) is mathematically interesting but in physics we do not encounter such
a force, except the $a=0$ limit.

From the expression of $V-E$ as given in Eq. (21) we get

$\displaystyle J(u)$ $\displaystyle\equiv$
$\displaystyle\frac{(V-E)}{L^{2}c^{2}}\frac{dV}{du}=\frac{d}{L^{2}c^{2}}(u+b)\,.$
(24)

yielding

$\displaystyle\frac{d^{2}u}{d\theta^{2}}+u=\frac{d}{L^{2}c^{2}}(u+b)\,,$ (25)

which gives the orbit equation of the relativistic particle which is acted on
by a force given by Eq. (23). As in general $d=L^{2}c^{2}(1-\zeta^{2})$ for a
stable closed orbit where $\zeta$ must be a rational number, we get

$\displaystyle\frac{1}{r}=\frac{1}{R}\cos(\zeta\theta)+\frac{b(1-\zeta^{2})}{\zeta^{2}}\,,$
(26)

where

$\displaystyle
R=Lc\zeta\left[\frac{b^{2}L^{2}c^{2}(1-\zeta^{2})}{\zeta^{2}}+a-m^{2}c^{4}\right]^{-1/2}\,.$
(27)

The equation of the orbit in Eq. (26) shows that in the most general case we
will have precession of the orbits dictated by the condition
$(2\pi+\delta\theta)\zeta=2\pi$, which predicts that the orbit precesses by an
angle

$\displaystyle\delta\theta=\frac{2\pi(1-\zeta)}{\zeta}\,,$ (28)

per orbit.

###  3.3 The case of large perturbations

Till now we have utilized first order perturbation from a circular orbit as
described in Eq. (10) in the beginning of this section. To include higher
order perturbations from a circular orbit we require more terms in the Taylor
series expansion of $J(u)$ in Eq. (10). The second order effects will come
from terms proportional to $(d^{2}J/du^{2})_{u_{0}}$. If the general form of
the potential $V(r)$ satisfies Eq. (21) then it is immediately clear from Eq.
(24) that all derivatives of $J(u)$, except the first, vanishes. Consequently
in the relativistic case it is impossible to restrict the constants $\zeta$,
$a$ and $b$ by higher order perturbation terms to the circular orbit. For
higher order perturbations from circular orbits the form the potential as
given in Eq. (21) remains the same.

##  4 Connection of the present work with some related works

One of the findings of the present article is related to the absence of
stable, circular orbits for relativistic particles in presence of a harmonic
oscillator potential. The trajectories of relativistic particles in a three
dimensional harmonic oscillator potential has been studied previously by L.
Homorodean in Ref. [8]. The method followed in the referred work is completely
equivalent to the one followed in the present work. It is interesting to note
that in Homorodean’s analysis the general shape of the orbit in the
relativistic case is not an ellipse, or a circle, but a rosette shaped curve.
In presence of the oscillator potential the angular momentum of the orbiting
particle with a specific energy has an upper bound. The trajectory of the
relativistic particle can only be a circle when it has the highest angular
momentum for a fixed energy. In Ref. [8] the author does not give any
information about the stability of the orbits. In the non-relativistic limit
the orbit of the particle can be circular. In the light of the findings in
Ref. [8] of Homorodean the prediction of the absence of a stable, circular
orbit in the oscillator potential is a sensible result.

Bertrand’s theorem in non-relativistic classical mechanics has inspired some
authors to propose a space-time (a metric to be precise) where any bounded
trajectory of a particle is periodic in nature. This kind of a space-time is
named as Bertrand space-time. The works of Perlick, Ballesteros, Enciso,
Herranz and Ragnisco, in Refs. [9, 10], try to generalize the results of the
classical Bertrand’s theorem on a flat 3-space to a curved 3-manifold. In Ref.
[9] the author found that a specific form of a space-time which is
asymptotically flat can support Keplerian orbits. The asymptotically flat
Bertrand space cannot support closed trajectories expected in an oscillator
type of potential. One of the findings of the present article predicts that
even in flat space relativistic effects forbid closed, stable trajectories of
particles in presence of an oscillator potential.

##  5 Conclusion

The outline of the article is based on the well known Bertrand’s theorem on
central potentials and orbits of particles as described in most of the
classical mechanics books. Like the non-relativistic case the relativistic
particle’s orbit around a potential source takes place in a plane where the
angular momentum and presumably the energy of the orbiting particle remains
constant. The main difference between the non-relativistic orbits and
relativistic orbits crops up in the orbit equation itself. Unlike the non-
relativistic case in the relativistic case the orbit equation depends upon the
total energy of the particle. The main aim of the article was to find out
possible forms of central potentials which can produce stable circular orbits
for relativistic particles. The stability condition for the orbits can be
transformed to a non-linear differential equation for the central potential.
It is seen that one of the solutions of the non-linear differential equation
for the central potential is just the normal Coulomb potential. But relativity
affects the properties of the orbits by curtailing the angular momentum values
beyond a certain limit. Except the Coulomb potential solution we find that the
stability equation has another general mathematically interesting solution
which is unlike any potential which we use in conventional physics. In a
specific limit the general solution reproduces the Newtonian or Coulomb form.
In the relativistic version of the central force problem we lack some
restrictions on the potential which can produce stable circular orbits. In the
non-relativistic version minimizing the effective potential one can figure out
the radius of the circular orbits and higher order perturbation corrections to
the stability condition of the orbits could be used for unravelling the exact
nature of the potential. In the relativistic version none of those
restrictions remain and consequently the general solution of the potential
contains some constants whose values cannot be analytically calculated.

An important fact which comes out from the article is about the non-existence
of the harmonic oscillator potential as a solution of the stability equation.
In non-relativistic treatment of the Bertrand’s theorem it is well known that
only two kinds of potentials can produce stable circular orbits, one is of the
Coulomb type and the other is of the harmonic oscillator type. The Coulomb
form of the potential passes the stability test for circular orbits but the
harmonic type does not.

The article was focussed on some mathematical properties of relativistic
particle orbits in a central potential. Before we finally conclude it is
pertinent to say some thing on the practical side of the relativistic central
force problem. Atomic physics always remains a store house of exciting
phenomena and one of the places where one may like to apply the tools of
relativistic central force problems lies inside the atom. This fact was
discussed in Ref. [2]. People have studied about the Schrödinger equation and
the Dirac equation in presence of the Coulomb potential. It can be quite
interesting to study the analogous problems using the potential presented in
this article instead of the Coulomb potential. This attempt can seriously shed
some light on the physics of the atoms.

More over as there exists some work on the general relativistic generalization
of Bertrand’s theorem one may expect that in the simplest of the situations,
where space-time remains flat, the results of the present work can be applied
for the orbits of very fast moving bodies interacting via Newtonian gravity
with a massive source. Acknowledgement: The authors acknowledge the
illuminating comments from H. C. Verma after he read the initial manuscript.
Many of his suggestions were implemented while preparing the final version of
the manuscript.

## References

  * [1] H. Goldstein, “Classical Mechanics”, $2^{\rm nd}$ edition, Narosa Publishing House, (1993)
  * [2] T. H. Boyer, Am. J. Phys. 72, 992-997, (2004)
  * [3] U. Torkelsson, Eur. J. Phys 19, 459-464, (1998)
  * [4] Z. Reut, [Q. J] Mech. Appl. Math., 39, 417-423, (1986)
  * [5] H. Frommert, Int. J. Theor. Phys., 35, 2631-2643, (1996)
  * [6] L. D. Landau, E. M. Lifshitz, “The classical theory of fields”, $4^{\rm th}$ edition, Elsevier, (2008)
  * [7] J. D. Jackson, “Classical Electrodynamics”, $2^{\rm nd}$ Edition, Wiley Eastern Limited, (1975).
  * [8] L. Homorodean, Europhys. Lett., 66, 8-13, (2004)
  * [9] V. Perlick, Class. Quantum Grav., 9, 1009-1021, (1992)
  * [10] A. Ballesteros, A. Enciso, F. J. Herranz and O. Ragnisco, Class. Quantum Grav., 9, 165005 (13pp) (2008)


# Testing trajectory-based determinism via time probability distributions

M. V. Scherer  A. D. Ribeiro  R. M. Angelo  Departamento de Física,
Universidade Federal do Paraná, 81531-990, Curitiba, PR, Brazil

###### Abstract

It is notorious that quantum mechanics (QM) cannot predict well-defined values
for all physical quantities. Less well-known, however, is the fact that QM is
unable to furnish probabilistic predictions even in emblematic scenarios such
as the double-slit experiment. In contrast, equipped with postulate
trajectories, Bohmian mechanics (BM) has inherited more predictive power. It
follows that, contrary to common belief, QM and BM are not just different
interpretations but distinct theories. This work formalizes the aforementioned
assertions and illustrates them through three case studies: (i) free particle,
(ii) free fall under a uniform gravitational field, and (iii) the double-slit
experiment. Specifically, we introduce a prescription for constructing a
flight-time probability distribution within generic trajectory-equipped
theories. We then apply our formalism to BM and derive probability
distributions that are unreachable by QM. Our results can, in principle, be
tested against real experiments, thereby assessing the validity of Bohmian
trajectories.

##  I Introduction

There are numerous physical and philosophical discussions related to quantum
mechanics (QM). However, one fact remains undeniable: the quantum formalism
has consistently succeeded in describing the distribution of detector clicks
in experiments. This success might lead some to believe that it provides
objective answers (though not interpretative) for all observable phenomena.
However, a critical examination of some traditional experiments reveals that
this is not always true. For example, consider a double-slit experiment [1, 2,
3, 4, 5, 6, 7], where the position $y$ of a particle is measured on a
detection screen located at $x=x_{\text{d}}$, with the sub-index “d” referring
to the detector. Typically, the particle’s flight-time is neither controlled
nor measured in such experiments. After many runs of the experiment, the
observer constructs the experimental time-independent probability density
${\wp_{y|x}}^{\mkern-22.0mu\text{exp}}(y|x_{\text{d}})$, which describes the
possible locations $y$ reachable by the particle, conditioned on the detection
line $x=x_{\text{d}}$. Now, how can one apply Born’s rule from the theoretical
$|\psi(x,y,t)|^{2}$ to fit the experimental data?

A first move towards answering this question consists of recognizing that
$|\psi(x,y,t)|^{2}$ is not a joint probability distribution for $x$, $y$, and
$t$. Instead, it refers to the probability density of finding the system at
the location $(x,y)$ given that time is guaranteed to be $t$. In this case,
the notation $\wp_{xy|t}(x,y|t)$111To shorten the notation, whenever the sub-
index of the symbol representing the probability density function repeats in
its argument, we may omit the latter. For example, $\wp_{xy|t}(x,y|t)$ would
be simply $\wp_{xy|t}$., emphasizing the conditioning on $t$, is preferred [8,
9, 10, 11, 12]. Rather than just allegorical rephrasing, this serves to stress
that the distribution has to be built for a single well-defined value of time.
Next, we can condition on $x=x_{\text{d}}$:

$\displaystyle\wp_{y|xt}(y|x_{\text{d}},t)=\frac{\wp_{xy|t}(x_{\text{d}},y|t)}{\wp_{x|t}(x_{\text{d}}|t)}=\frac{|\psi(x_{\text{d}},y,t)|^{2}}{\int\\!dy\,|\psi(x_{\text{d}},y,t)|^{2}}.$
(1)

Unless one post-selects the experimental data to a particular value of $t$,
which actually is never done in practice, we still do not have a formal
comparison with ${\wp_{y|x}}^{\mkern-22.0mu\text{exp}}(y|x_{\text{d}})$. Of
course, our first impulse would be to assign some mean value for time in Eq.
(1) or even to integrate over some time window, but these would be extra
assumptions. This underscores that experimental results cannot always be
straightforwardly confronted with QM’s predictions.

Experiments like the ones quoted above are usually carried out in the far-
field regime, so semiclassical analysis generally applies [13, 14, 15, 16].
Nevertheless, recent advancements in ultra-fast detector technology [17, 18,
19, 20] have made it possible to explore regimes where semiclassics is
insufficient to describe the system’s behavior [13, 21]. Therefore, it is
mandatory to develop a strategy to deal with time. A viable route to this end
would be to search for a probability distribution for the flight-time,
$\wp_{t|x}(t|x_{\text{d}})$, and then use probability theory to obtain the
desired distribution,

$\displaystyle\wp_{y|x}(y|x_{\text{d}})=\int\\!dt\,\wp_{y|xt}(y|x_{\text{d}},t)\,\wp_{t|x}(t|x_{\text{d}}).$
(2)

It turns out, however, that the formalism of QM does not provide any
prescription for the construction of $\wp_{t|x}(t|x_{\text{d}})$, thus
remaining silent regarding the time distribution of clicks along the detection
line $x=x_{\text{d}}$.

As a first contribution of this article, we demonstrate that trajectory-
equipped theories can naturally address this issue. After presenting a general
argument, our attention will be directed towards Bohmian mechanics (BM) [22,
23], a model that reinstates determinism while upholding probabilistic
predictions consistent with QM. Although different from ours, approaches are
already known that use Bohmian trajectories to construct time probability
distributions [24, 25, 26]. In fact, several trajectory-based models exist
[27, 28, 29] that allow access to the probability distribution of time [30,
31, 32]. Other approaches to time in QM are worth citing. While the Kijowski
method [33, 34] employs semiclassical approximations to derive temporal
predictions, the Page-Wooters mechanism [35] considers time as an inaccessible
coordinate that arises from correlations between subsystems of the global
state. Similar approaches conceiving time as an emergent property of physical
systems can be found in Refs. [36, 37, 38]. The spacetime-symmetric extension
of QM [9] treats time on par with space, resulting in a Schrödinger equation
for time. This approach has found application to the problem of the tunneling
time [39]. Finally, there are phenomenological modeling approximations of the
detection process that yield probability distributions of time, such as the
absorbing boundary rule [40, 41, 42, 43], path-integral with absorbing
boundary [44, 45, 46], and others [47, 48, 49, 50, 51]. Each of these
approximations and formalisms can be tested using equation (2) as they access
the probability distributions $\wp_{y|xt}$ and $\wp_{t|x}$, with the sole
requirement being experimental data for comparison, without the need for
experiments involving direct time measurements.

In this work, a formalism based on the hypothesis of trajectory-based
determinism will be developed to derive spatial and temporal probability
distributions. Supplemented with minimal assumptions, this formalism enables
the computation of probability distributions $\wp_{x|t}$ and $\wp_{t|x}$, and
finally the one in Eq. (2). We then employ BM to give a concrete example of
the formalism usage and compare our results with those of other approaches
[24, 25, 52, 53]. Besides this introduction, the text is organized as follows.
In Sec. II, we introduce the basis for our model, which involves two
postulates and the development of conditional probability distributions, one
of them for the flight-time. In Sec. III, we apply our formalism to some
emblematic physical problems and compare our results with other approaches
(when available). We close the article in Sec. IV with our final remarks.

##  II Trajectory-equipped theory (TET)

In this section, we present our formalism, treating separately two types of
deterministic trajectories, namely, those defined in terms of a velocity field
and those in phase space. To simplify the discussion, we restrict ourselves to
the one-dimensional case.

###  II.1 Trajectories associated with a velocity field

In this scenario, trajectories emerge from a first-order differential equation
$\dot{q}=f(q,t)$, where $f$ is some real function representing the velocity
field. Setting $q_{0}$, the coordinate of the particle at $t=0$, the
coordinate $q$ of the particle at the instant $t$ is uniquely defined. As is
the case in BM, $q_{0}$ is the only boundary condition needed, and the
mechanical momentum of a particle of mass $m$ is just given by $p=mf(q,t)$.
Now, as far as the triple $\\{q,q_{0},t\\}$ is concerned, here is the catch:
just as the specification of $q_{0}$ and $t$ implies $q$, the specification of
$q$ and $q_{0}$ (respectively, $q$ and $t$) should imply $t$ (respectively,
$q_{0}$). In addition, due to the functional structure connecting the triple,
one variable inherits the uncertainties of the others, provided physical
principles are not violated. We rigorously formalize these ideas as follows:

  * •

Determinism hypothesis.—There exists a smooth trajectory $q=Q(q_{0},t)$
connecting an initial coordinate $q_{0}$ to a point $q$ over a time lapse $t$.
For bounded dynamics, $q$ is expected to be a recurrence point, reachable at
several instants $t_{j}\in\tau_{n}\equiv\\{t_{1},t_{2},\dots,t_{n}\\}$.
Additionally, in typical well-behaved dynamics, $Q$ is invertible (at least in
parts). Therefore, there must be functions $Q_{0}$ and
$\\{T_{\mkern-4.0muj}\\}_{j=1}^{n}$ allowing us to express determinism in
terms of the relations

$\displaystyle q=Q(q_{0},t),\quad q_{0}=Q_{0}(q,t),\quad
t_{j}=T_{\mkern-4.0muj}(q,q_{0}),$ (3)

with $t\in\tau_{n}$. In this trajectory-equipped theory (TET), probability
distributions are given by

$\displaystyle\wp_{q|q_{0}t}(q|q_{0},t)=\delta\big{(}q-Q(q_{0},t)\big{)},$
(4a)
$\displaystyle\wp_{t|qq_{0}}(t|q,q_{0})=\frac{1}{n}\sum_{j=1}^{n}\delta\big{(}t-T_{\mkern-4.0muj}(q,q_{0})\big{)},$
(4b)

where $\delta$ denotes the Dirac delta function. As expected, here we have
non-negative distributions satisfying
$\int\\!dq\,\wp_{q|q_{0}t}(q|q_{0},t_{j})=\int_{0}^{\infty}\\!dt\,\wp_{t|qq_{0}}(t|q,q_{0})=1$.

  * •

Preparation independence hypothesis.— The probability distribution for $q_{0}$
derives from some experimental procedure that is in no way physically
conditioned to $q$ and $t>0$, hence

$\displaystyle\wp_{q_{0}|qt}(q_{0}|q,t)=\wp_{q_{0}}(q_{0})\equiv\rho_{0}(q_{0}).$
(5)

The preparation $\rho_{0}(q_{0})$ should be specified in each statistical
theory. In QM, $\rho_{0}(q_{0})=\langle q_{0}|\hat{\rho}_{0}|q_{0}\rangle$,
with $\hat{\rho}_{0}$ the density operator at $t=0$.

In what follows, we will build upon the previously mentioned assumptions,
using the definition of marginal probability,

$\displaystyle\wp_{v|r}(v|r)=\int\\!du\,\wp_{v|ru}(v|r,u)\,\wp_{u|r}(u|r),$
(6)

along with the property of the Dirac delta function,

$\displaystyle\delta\left(g(u)\right)=\sum_{i=1}^{N}\frac{\delta(u-u_{i})}{|\partial_{u}g(u)|_{u=u_{i}}},$
(7)

where $\partial_{u}g\equiv\partial g/\partial u$ and $g(u_{i})=0$ for $1\leq
i\leq N$. Plugging the distributions (4a) and (5) into Eq. (6), and solving
the integral over $q_{0}$ with the property (7), we obtain

$\displaystyle\wp_{q|t}(q|t)=\left[\frac{\rho_{0}(q_{0})}{\left|\partial_{q_{0}}Q(q_{0},t)\right|}\right]_{q_{0}=Q_{0}(q,t)}.$
(8)

The modulus operation appearing in the denominator guarantees that
$\wp_{q|t}\geq 0$. The probability distribution for flight-times is
constructed similarly. Using the relations (4b)-(7), we end up with

$\displaystyle\wp_{t|q}(t|q)=\frac{1}{n}\sum_{j=1}^{n}\left[\frac{\rho_{0}\big{(}q_{0}\big{)}}{\left|\partial_{q_{0}}T_{\mkern-4.0muj}(q,q_{0})\right|}\right]_{q_{0}=Q_{0j}(q,t)},$
(9)

where the function $Q_{0j}(q,t)$ is achieved from manipulating the equation
$t=T_{\mkern-4.0muj}(q,q_{0})$ so that $\smash{q_{0}=Q_{0j}(q,t)}$.

For future reference, we now indicate how the probability distributions of
flight-times and positions can be associated with each other, in particular in
the BM context. Starting with $\wp_{q|q_{0}t}$, we can use the property (7) to
write

$\displaystyle\wp_{q|q_{0}t}=\delta\big{(}q-Q(q_{0},t)\big{)}=\sum_{j=1}^{n}\frac{\delta\big{(}t-T_{\mkern-4.0muj}(q,q_{0})\big{)}}{v_{j}(q,q_{0})},$
(10)

where we have introduced
$v_{j}(q,q_{0})\equiv|\partial_{t}Q(q_{0},t)|_{t=T_{\mkern-4.0muj}(q,q_{0})}\geq
0$. Multiplying Eq. (10) by $\rho_{0}(q_{0})$ and integrating over $q_{0}$
gives

$\displaystyle\wp_{q|t}(q|t)=\sum_{j=1}^{n}\left[\frac{\rho_{0}(q_{0})}{v_{j}(q,q_{0})\left|\partial_{q_{0}}T_{\mkern-4.0muj}(q,q_{0})\right|}\right]_{q_{0}=Q_{0j}(q,t)}.$
(11)

In the specific case where $v_{j}(q,q_{0})\equiv v(q,t)$, independently of
$j$, we find, from Eq. (9),

$\displaystyle\wp_{t|q}(t|q)=\frac{v(q,t)\,\wp_{q|t}(q|t)}{n}.$ (12)

Specializing to BM, $v(q,t)$ is readily identified with the velocity field,
and $|J(q,t)|\equiv v(q,t)\,\wp_{q|t}(q|t)$ with the absolute value of the
probability current density, so we have

$\displaystyle{\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}}(t|q)=\frac{|J(q,t)|}{n}.$
(13)

This formula is a generalization of the typical one [24, 25], valid only when
$q$ is not a recurrence point $(n=1)$. In fact, it is crucial to realize that
our approach generalizes the existing one in two aspects. First, it applies to
all trajectory-equipped theories, not only to BM. Second, it naturally
accommodates bounded oscillatory dynamics, an issue that, to our knowledge,
has remained unaddressed so far.

###  II.2 Phase space trajectories

We now consider trajectories defined in the phase space $qp$ through real
functions $q=\mathcal{Q}(\mathbf{r}_{0},t)$ and
$p=\mathcal{P}(\mathbf{r}_{0},t)$, constrained by the initial condition
$\mathbf{r}_{0}=(q_{0},p_{0})$ at $t=0$. Importantly, since $q$ and $p$ are
treated as independent canonical variables, the uncertainties associated with
$q_{0}$ and $p_{0}$ can similarly be independent. The initial probability
distribution is denoted as $\rho_{0}(\mathbf{r}_{0})$ or, equivalently,
$\rho_{0}\big{(}q_{0},p_{0}\big{)}$. For simplicity, in what follows, we
assume that there are no recurrences (i.e., the trajectory passes through the
point $\mathbf{r}=(q,p)$ only once), so the inversion of the trajectories
yields only the functions $t=\mathcal{T}_{\mkern-4.0muq}(q,\mathbf{r}_{0})$
and $t=\mathcal{T}_{\mkern-5.0mup}(p,\mathbf{r}_{0})$. Again, we propose an
approach relying on two fundamental assumptions.

  * •

Determinism hypothesis.—There exists a trajectory in phase space defined by
smooth functions $q=\mathcal{Q}(\mathbf{r}_{0},t)$ and
$p=\mathcal{P}(\mathbf{r}_{0},t)$, connecting an initial point
$\mathbf{r}_{0}=(q_{0},p_{0})$ to a final point $\mathbf{r}=(q,p)$ over a time
lapse $t$. Assuming that the involved functions are invertible and that there
are no recurring points, the determinism hypothesis can be expressed by the
relations

$\displaystyle\begin{array}[]{lll}q=\mathcal{Q}(\mathbf{r}_{0},t),&&q_{0}=\mathcal{Q}_{0q}(q,p_{0},t),\\\\[4.0pt]
p_{0}=\mathcal{P}_{0q}(q,q_{0},t),&&t=\mathcal{T}_{\mkern-4.0muq}(q,\mathbf{r}_{0}),\end{array}$
(16)

and

$\displaystyle\begin{array}[]{lll}p=\mathcal{P}(\mathbf{r}_{0},t),&&q_{0}=\mathcal{Q}_{0p}(p,p_{0},t),\\\\[4.0pt]
p_{0}=\mathcal{P}_{0p}(p,q_{0},t),&&t=\mathcal{T}_{\mkern-5.0mup}(p,\mathbf{r}_{0}).\end{array}$
(19)

Note that the symbol $\mathcal{Q}_{0q}$ (respectively, $\mathcal{P}_{0q}$)
refers to a function originating from the inversion of
$q=\mathcal{Q}(\mathbf{r}_{0},t)$, with output $q_{0}$ (respectively,
$p_{0}$). Similar reasoning has motivated us to name the other functions
appearing in Eqs. (16) and (19). In addition, note that the pairs
$\big{\\{}\mathcal{T}_{\mkern-4.0muq},\mathcal{T}_{\mkern-5.0mup}\big{\\}}$,
$\big{\\{}\mathcal{Q}_{0q},\mathcal{Q}_{0p}\big{\\}}$, and
$\big{\\{}\mathcal{P}_{0q},\mathcal{P}_{0p}\big{\\}}$ consist of distinct
functions. Based on these relations, we propose the following conditional
probabilities:

$\displaystyle\wp_{qp|\mathbf{r}_{0}t}=\delta\big{(}q-\mathcal{Q}(\mathbf{r}_{0},t)\big{)}\,\delta\big{(}p-\mathcal{P}(\mathbf{r}_{0},t)\big{)},$
(20a)
$\displaystyle\wp_{qt|p\mathbf{r}_{0}}=\delta\big{(}q-\mathcal{Q}(\mathbf{r}_{0},t)\big{)}\,\delta\big{(}t-\mathcal{T}_{\mkern-5.0mup}(p,\mathbf{r}_{0})\big{)},$
(20b)
$\displaystyle\wp_{pt|q\mathbf{r}_{0}}=\delta\big{(}p-\mathcal{P}(\mathbf{r}_{0},t)\big{)}\,\delta\big{(}t-\mathcal{T}_{\mkern-4.0muq}(q,\mathbf{r}_{0})\big{)}.$
(20c)

The approach is such that starting with two postulated functions $\mathcal{Q}$
and $\mathcal{P}$ involving the variables $\\{q,p,q_{0},p_{0},t\\}$, one can
construct a joint probability distribution for two of them out of the
specification of the others. Note that while the distribution (20a) admits the
interpretation
$\wp_{qp|\mathbf{r}_{0}t}=\wp_{q|\mathbf{r}_{0}t}\,\wp_{p|\mathbf{r}_{0}t}$,
meaning that $q$ and $p$ are uncorrelated variables both conditioned to
$\\{\mathbf{r}_{0},t\\}$, the same does not hold for the ones involving time.
Equation (20b), for instance, can be interpreted as
$\wp_{qt|p\mathbf{r}_{0}}=\wp_{q|\mathbf{r}_{0}t}\,\wp_{t|p\mathbf{r}_{0}}$,
with the probability for one of the variables, $\wp_{q|\mathbf{r}_{0}t}$,
depending explicitly on the other. All this is consistent with the premise
that $q$ and $p$ are independent canonical variables, while $q$ and $t$ are
not.

  * •

Preparation independence hypothesis.—The joint probability for
$\mathbf{r}_{0}$ derives from some experimental procedure that is in no way
physically conditioned to $\mathbf{r}$ and $t>0$, hence

$\displaystyle\wp_{\mathbf{r}_{0}|\mathbf{r}\mkern
1.0mut}(\mathbf{r}_{0}|\mathbf{r},t)=\wp_{\mathbf{r}_{0}}(\mathbf{r}_{0})\equiv\rho_{0}(\mathbf{r}_{0}).$
(21)

The preparation $\rho_{0}(\mathbf{r}_{0})$ should be specified in each
statistical theory. In QM, one possibility is to use the Husimi distribution
[54], $\langle\alpha|\hat{\rho}_{0}|\alpha\rangle/\pi$, where $\ket{\alpha}$
denotes a coherent state with complex label $\alpha=\alpha(\mathbf{r}_{0})$
parametrized in terms of phase-space coordinates.

Utilizing the same procedure as before, we integrate the product of
distributions (20a) and (21) over $\mathbf{r}_{0}$ to obtain, with the aid of
property (7),

$\displaystyle\wp_{qp|t}(q,p|t)=\left[\frac{\rho_{0}\big{(}\mathcal{Q}_{0q}(q,p_{0},t),p_{0}\big{)}}{\mathcal{D}_{q_{0}}^{\mathcal{Q}}(q,p_{0},t)~{}\mathcal{D}_{p_{0}}^{\mathcal{P}}(q,p_{0},t)}\right]_{p_{0}=\mathcal{P}_{\mkern-6.0mu\mathcal{Q}}(q,p,t)},$
(22)

where

$\displaystyle\begin{aligned}
\mathcal{D}_{q_{0}}^{\mathcal{Q}}(q,p_{0},t)&\equiv\left|\partial_{q_{0}}\mathcal{Q}\big{(}q_{0},p_{0},t\big{)}\right|_{q_{0}=\mathcal{Q}_{0q}(q,p_{0},t)},\\\
\mathcal{D}_{p_{0}}^{\mathcal{P}}(q,p_{0},t)&\equiv\left|\partial_{p_{0}}\mathcal{P}\big{(}\mathcal{Q}_{0q}(q,p_{0},t),p_{0},t\big{)}\right|,\end{aligned}$
(23)

and $p_{0}=\mathcal{P}_{\mkern-6.0mu\mathcal{Q}}(q,p,t)$ is found by isolating
$p_{0}$ in the equation
$p=\mathcal{P}\big{(}\mathcal{Q}_{0q}(q,p_{0},t),p_{0},t\big{)}$. Via similar
procedures, from Eq. (20b), we obtain

$\displaystyle\wp_{qt|p}(q,t|p)=\left[\frac{\rho_{0}\big{(}\mathcal{Q}_{0q}(q,p_{0},t),p_{0}\big{)}}{\mathcal{D}_{q_{0}}^{\mathcal{Q}}(q,p_{0},t)~{}\mathcal{D}_{p_{0}}^{\mathcal{T}}(q,p,p_{0},t)}\right]_{p_{0}=\mathcal{P}_{\mathcal{T}}(q,p,t)},$
(24)

where

$\displaystyle\mathcal{D}_{p_{0}}^{\mathcal{T}}(q,p,p_{0},t)\equiv\left|\partial_{p_{0}}\mathcal{T}_{\mkern-5.0mup}\big{(}p,\mathcal{Q}_{0q}(q,p_{0},t),p_{0}\big{)}\right|$
(25)

and $p_{0}=\mathcal{P}_{\mathcal{T}}(q,p,t)$ is found by isolating $p_{0}$ in
the equation
$t=\mathcal{T}_{\mkern-5.0mup}\big{(}p,\mathcal{Q}_{0q}(q,p_{0},t),p_{0}\big{)}$,
and, from Eq. (20c),

$\displaystyle\wp_{pt|q}(p,t|q)=\left[\frac{\rho_{0}\big{(}\mathcal{Q}_{0q}(q,p_{0},t),p_{0}\big{)}}{\mathcal{D}_{q_{0}}^{\mathcal{T}}(q,p_{0},t)~{}\mathcal{D}_{p_{0}}^{\mathcal{P}}(q,p_{0},t)}\right]_{p_{0}=\mathcal{P}_{\mkern-6.0mu\mathcal{Q}}(q,p,t)},$
(26)

where function $\mathcal{P}_{\mkern-6.0mu\mathcal{Q}}$ is same one that
appears in Eq. (22) and

$\displaystyle\mathcal{D}_{q_{0}}^{\mathcal{T}}(q,p_{0},t)\equiv\left|\partial_{q_{0}}\mathcal{T}_{\mkern-4.0muq}\big{(}q,q_{0},p_{0}\big{)}\right|_{q_{0}=\mathcal{Q}_{0q}(q,p_{0},t)}.$
(27)

The equations (22), (24), and (26) provide instructions for deriving
probability distributions for a trajectory in phase space. To find them,
specifically in the form of $\wp_{q|t}$ and $\wp_{t|q}$, marginalization
should be performed with respect to the momentum $p$.

##  III Case studies

In this section, we apply our formalism to simple physical systems and, when
pertinent, we compare our results with those of other approaches.

###  III.1 Free particle

We start considering the one-dimensional dynamics of a free particle with mass
$m$, prepared in a Gaussian state written as
$\psi(q,t=0)=\sqrt{G_{q}(\bar{q}_{0},\sigma_{0})}\,\mbox{e}^{i\bar{p}_{0}\left(q-\bar{q}_{0}\right)/\hbar}$,
where

$\displaystyle
G_{u}(\bar{u},\Delta)\coloneqq\frac{1}{\sqrt{2\pi\Delta^{2}}}\exp{\left[-\frac{(u-\bar{u})^{2}}{2\Delta^{2}}\right]}.$
(28)

The parameter $\bar{q}_{0}$ ($\bar{p}_{0}$) is the mean position (momentum)
and $\sigma_{0}$ is the spatial uncertainty. For the role of TET, we choose
two participants: classical statistical mechanics (or Liouvillian mechanics)
and BM. For the classical model, the phase-space trajectory is given by

$\displaystyle q$
$\displaystyle=\mathcal{Q}(\mathbf{r}_{0},t)=q_{0}+\frac{p_{0}t}{m},$
$\displaystyle p$ $\displaystyle=\mathcal{P}(\mathbf{r}_{0},t)=p_{0},$

while
$\rho_{0}(\mathbf{r}_{0})=G_{q_{0}}\left(\bar{q}_{0},\sigma_{0}\right)\,G_{p_{0}}\left(\bar{p}_{0},\frac{\hbar}{2\sigma_{0}}\right)$
stands for the phase-space initial distribution. The joint probability
distribution is then obtained using equation (22). The result reads

$\displaystyle{\wp_{qp|t}}^{\scriptscriptstyle\mkern-40.0mu\text{C}\mkern
28.0mu}(q,p|t)=\rho_{0}\bigg{(}q-{\textstyle\frac{pt}{m}},p\bigg{)}=G_{q}\bigg{(}\bar{q}_{t},\sigma_{0}\bigg{)}\,G_{p}\bigg{(}\bar{p}_{0},{\textstyle\frac{\hbar}{2\sigma_{0}}}\bigg{)},$

where $\bar{q}_{t}\equiv\bar{q}_{0}+\frac{pt}{m}$ and “$\mathrm{C}$” stands
for “classical”. Marginalizing over $p$ gives

$\displaystyle{\wp_{q|t}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern
13.0mu}(q|t)=\int\\!dp\,{\wp_{qp|t}}^{\scriptscriptstyle\mkern-40.0mu\text{C}\mkern
28.0mu}(q,p|t)=G_{q}\left(\bar{q}_{0}+\frac{\bar{p}_{0}t}{m},\sigma_{t}\right),$
(29)

where

$\displaystyle\sigma_{t}=\sigma_{0}\sqrt{1+(\omega
t)^{2}},\qquad\omega=\hbar/(2m\sigma_{0}^{2}).$ (30)

As expected, these results plainly match the ones obtained from the
Liouvillian formalism, where
${\wp_{qp|t}}^{\scriptscriptstyle\mkern-40.0mu\text{C}\mkern 28.0mu}$ emerges
from the application of the inverse map of the trajectory to the initial
probability distribution [55]. In addition,
${\wp_{q|t}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern
13.0mu}=|\psi(q,t)|^{2}$, meaning that the classical distribution mirrors the
quantum one, as already known for the free particle system.

To compute the flight-time probability distribution, we start from Eq. (26) to
obtain

$\displaystyle{\wp_{pt|q}}^{\scriptscriptstyle\mkern-40.0mu\text{C}\mkern
28.0mu}(p,t|q)=\frac{|p|}{m}\,G_{q}\bigg{(}\bar{q}_{t},\sigma_{0}\bigg{)}\,G_{p}\bigg{(}\bar{p}_{0},{\textstyle\frac{\hbar}{2\sigma_{0}}}\bigg{)}.$
(31)

Now, to get ${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern
13.0mu}$, one needs to integrate this result over $p$. Due to its complexity,
the analytical result will be omitted, but it is numerically illustrated in
Fig. 1, constrained to the detection position $q=0$. This graph will be
discussed later when confronted with the next approaches.

Figure 1: Probability distributions
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern 13.0mu}(t|0)$
(thick yellow line),
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern 1.0mu}(t|0)$
(thin red line), and
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{K}\mkern 13.0mu}(t|0)$
(dashed black line), for the flight-time of a free particle detected at the
position $q=0$, according to two trajectory-equipped theories, the classical
one (C) and Bohmian mechanics (BM), and the Kijowski method (K). Panels (a)
and (b) depict distributions for two parameter sets, namely,
${\\{\sigma_{0}=2,\bar{q}_{0}=-15,\bar{p}_{0}=1,m=0.5,\hbar=1\\}}$ and
${\\{\sigma_{0}=2,\bar{q}_{0}=-15,\bar{p}_{0}=20,m=10,\hbar=1\\}}$,
respectively. All physical quantities and distributions are in arbitrary
units. The gray vertical lines denote the flight-time
$-m\bar{q}_{0}/\bar{p}_{0}$ for the center of the classical distribution,
while the dashed green lines represent the average times, which converge to
the same value for all distributions. A complete agreement is observed among
the models.

As a second TET, we consider BM. According to Bohmian theory, the trajectory
for the free particle is given by [24]

$\displaystyle Q(q_{0},t)=\bar{q}_{0}+\frac{\bar{p}_{0}t}{m}+\sqrt{1+(\omega
t)^{2}}\left(q_{0}-\bar{q}_{0}\right).$ (32)

With this function and its inverse $Q_{0}(q,t)$, we come to Eq. (8) to show
that ${\wp_{q|t}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern
1.0mu}={\wp_{q|t}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern
13.0mu}=|\psi(q,t)|^{2}$. Through Eq. (9), we derive the flight-time
distribution:

$\displaystyle{\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern
1.0mu}(t|q)=\left|\frac{\bar{p}_{0}}{m}+\frac{(\omega t)^{2}}{1+(\omega
t)^{2}}\left(\frac{q-\bar{q}_{0}}{t}-\frac{\bar{p}_{0}}{m}\right)\right|\,|\psi(q,t)|^{2},$
(33)

where the term inside the modulus is the velocity field $v(q,t)$.

We also bring to the discussion the distribution deriving from the Kijowski
method [33, 34], which relies on the hypothesis of a generalized measurement
operator for flight-time. Specifically for the arrival location $q=0$ and
Gaussian state, the resulting distribution reads [34]

where $\Theta$ is the step function.

In Fig. 1, we depict the distributions
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{C}\mkern 13.0mu}(t|0)$,
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern 1.0mu}(t|0)$, and
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{K}\mkern 13.0mu}(t|0)$ for
two sets of parameters. The first remark is that all theories yield precisely
the same distribution. To a certain extent, this is not surprising, as the
model under inspection is the simplest possible: the dynamics is linear, and
the preparation is Gaussian. On the other hand, there is subtler information
coming from the results. In the two panels considered, the mean velocities are
identical, suggesting that the mean flight-times should be equal as well.
However, this expectation is not confirmed due to the disparity in the masses.
In panel (a), the mass is significantly smaller than the mass in panel (b),
resulting in a larger spreading factor $\omega t$ and, as a consequence, an
enhancement of the temporal asymmetry induced the term $(\omega
t)^{2}/[1+(\omega t)^{2}]$. This effect can be clearly appreciated in the BM
formalism through the relations (32) and (33). As illustrated in panel (b),
the asymmetry in the time-probability distribution disappears as $\omega\to
0$.

###  III.2 Free fall under a uniform gravitational field

Let’s delve into a scenario involving a freely falling particle, where the
trajectory intersects the vertical coordinate $q=0$ at two distinct instances.
As discussed previously, the TET posits that the probability distribution for
flight-times is governed by Eq. (9) or, upon marginalization over $p$, Eq.
(26). In the present subsection, however, our study will be restricted to
trajectories generated by a velocity field, implying that the latter equation
does not apply. In addition, we will set TET to be BM. With this choice, we
remind that the Muga-Leavens (ML) result is retrieved [24], with the caveat
that it falls short in describing some systems involving recurrences (the
back-flow problem). Here, the idea is to show that our formalism performs
better than ML’s in characterizing the flight-time probability distribution
when this behavior manifests. Given our purposes, the analysis of the free
fall motion will be effectively restricted to one dimension, where we can
directly apply our formalism.

Figure 2: Bohmian trajectories of a freely falling particle for an ensemble of
100 initial positions $q_{0}$. The number of initial conditions located inside
the interval $[q_{0},q_{0}+dq_{0}]$ is proportional to
$|\psi(q_{0},0)|^{2}dq_{0}$, to be statistically representative. Panels (a)
and (b) depict the distributions for two sets of parameters, namely,
${\\{\sigma_{0}=1,\bar{q}_{0}=-8,\bar{p}_{0}=10,m=0.5,g=10,\hbar=1\\}}$ and
${\\{\sigma_{0}=2,\bar{q}_{0}=-8,\bar{p}_{0}=40,m=2,g=10,\hbar=1\\}}$,
respectively. The gray line $q=0$ indicates the position where the detector is
located. All physical quantities are given in arbitrary units.

Consider a particle of mass $m$ immersed in a uniform gravitational field and
prepared in a Gaussian state expressed by
$\psi(q,t=0)=\sqrt{G_{q}(\bar{q}_{0},\sigma_{0})}\,\mbox{e}^{i\bar{p}_{0}(q-\bar{q}_{0})/\hbar}$,
where $\bar{q}_{0}$ ($\bar{p}_{0}$) is the mean position (momentum) and
$\sigma_{0}$ is the initial spatial uncertainty. The Bohmian trajectory for
this case is directly calculated as

$\displaystyle
Q(q_{0},t)=\bar{q}_{0}+\frac{\bar{p}_{0}t}{m}-\frac{gt^{2}}{2}+\sqrt{1+(\omega
t)^{2}}\,\left(q_{0}-\bar{q}_{0}\right),$ (34)

where $\omega$ is given in Eq. (30) and $g$ is the gravitational constant. In
what follows, the parameters will be chosen to ensure that the particle
crosses the coordinate $q=0$ at two distinct instants for practically all
significant initial conditions $q_{0}$. To illustrate this kind of dynamics,
in Fig. 2, an ensemble of Bohmian trajectories (34) is showcased for two
distinct sets of parameters. The choice of $q_{0}$ for each trajectory was
taken considering the initial distribution $\psi(q,0)$, namely, the density of
points $q_{0}$ is proportional to $|\psi(q_{0},0)|^{2}$. As demanded, all
trajectories shown in the graph have two passages through $q=0$.

The flight-time probability distribution for the ML formalism derives from the
probability current density [25]. Direct calculations yield

$\displaystyle{\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{ML}\mkern
1.0mu}(t|q)=\frac{1}{\mathcal{N}}|v(q,t)|\,|\psi(q,t)|^{2},$ (35)

where the velocity field takes the form

$\displaystyle v(q,t)=\frac{\bar{p}_{0}}{m}-gt+\frac{(\omega t)^{2}}{1+(\omega
t)^{2}}\left(\frac{q-\bar{q}_{0}}{t}-\frac{\bar{p}_{0}}{m}+\frac{gt}{2}\right),$
(36)

$\mathcal{N}$ is a normalization constant, and $\psi(q,t)$ is the solution to
the Schrödinger equation.

To apply our method (9), we first need to invert the function (34) to find
solutions $t_{j}=T_{\mkern-4.0muj}(q,q_{0})$, for $j=1,2$, corresponding to
the two instants at which the trajectory crosses $q=0$ (see Fig. 2). This
entails solving a quartic equation for time and then discarding two
nonphysical roots. To discern the physical solutions, we require that
$t_{j}=T_{\mkern-4.0muj}(q,q_{0})\geq 0$, for any $\\{q,q_{0}\\}$, and the
consistency condition $t=T_{\mkern-4.0muj}\big{(}q,Q_{0}(q,t)\big{)}$. Put
simply, utilizing the inversion of the initial position $q_{0}=Q_{0}(q,t)$,
which is unique given $q$ and $t$, allows for finding $q_{0}$. If this initial
condition is applied to the function $T_{\mkern-4.0muj}$ and yields $t_{j}=t$,
the solution is deemed physical for a given time instant $t$. This methodology
aids in identifying the physical solutions for each time interval. Although
these criteria suffice for finding the physical solutions in the current
scenario, in more intricate cases, additional criteria may be necessary. For
practical purposes, numerical calculations will be conducted to obtain the
inverse trajectories and their derivatives, implementing them in Eq. (9).

Figure 3: Flight-time probability distributions for the problem of a freely
falling particle, obtained through our formalism equipped with Bohmian
trajectories [${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern
1.0mu}(t|0)$, red solid line], the ML method
[${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{ML}\mkern 1.0mu}(t|0)$,
dashed-black line], and via direct counting of the number of Bohmian
trajectories (blue dots). For the latter, we time evolved $10^{4}$ initial
conditions distributed according to $|\psi(q_{0},0)|^{2}$, counting those
arriving at $q=0$ after a time lapse $t$. Panels (a) and (b) display
distributions for two sets of parameters, namely,
${\\{\sigma_{0}=1,\bar{q}_{0}=-8,\bar{p}_{0}=9,m=0.5,g=10,\hbar=1\\}}$ and
${\\{\sigma_{0}=2,\bar{q}_{0}=-8,\bar{p}_{0}=36,m=2,g=10,\hbar=1\\}}$,
respectively. All physical quantities are in arbitrary units.

In Fig. 3, the flight-time probability distribution is depicted using the ML
method (black dashed curve) and our formalism equipped with Bohmian
trajectories (red solid curve). These results are also compared with a direct
counting of the number of Bohmian trajectories (blue dots), according to the
following procedure. We constructed an ensemble of $10^{4}$ initial positions,
respecting the probability density of initial positions,
$|\psi(q_{0},0)|^{2}$. Then, we evolved each component of the ensemble using
the Bohmian trajectory (34). Each time a trajectory intersected $q=0$, we
recorded its flight-time, and, at the end, we built the dotted distribution
shown in Fig. 3. In panel (a), the agreement (disagreement) between our
formalism (ML approach) and the direct counting is notable. This shows that
the TET formalism is somewhat superior to the ML approach in describing
systems exhibiting recurrences. The peaks of
${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern 1.0mu}(t|0)$ are
located precisely over the most probable classical times, that is, the
instants for which the center of the Gaussian distribution reaches $q=0$. The
distribution ${\wp_{t|q}}^{\scriptscriptstyle\mkern-25.0mu\text{ML}\mkern
1.0mu}(t|0)$, on the other hand, fails to predict the second peak, the one for
the recurrence time. As we approach a more classical regime (large mass), the
velocity fields (36) at $t_{1,2}$ become more symmetric, in the sense that
$|v(q=0,t_{1})|=|v(q=0,t_{2})|$. In this regime, as discussed at the end of
subsection II.1, the ML approach becomes as effective as ours [see panel (b)].

###  III.3 Double-slit experiment

Referring back to the problem discussed in the introduction, we now apply our
formalism to provide a theoretical counterpart of
${\wp_{y|x}}^{\mkern-22.0mu\text{exp}}(y|x_{\text{d}})$ for the double-slit
experiment with particles, where the detection screen is situated in the line
$(x_{\text{d}},y)$. The slits occupy the positions $\pm{y}_{\text{s}}$ over
the line $(0,y)$. Here, we also use only BM as TET. Because the components $x$
and $y$ of the particle’s motion are decoupled from each other, we have at
hand two independent one-dimensional problems, which are connected only by the
variable time. For the $x$ direction, we take the usual wave-plane solution
with box normalization:

$\displaystyle\psi_{x}(x,t)=\begin{cases}x_{\text{d}}^{-1/2}\mbox{e}^{i(\bar{p}_{0x}x-Et)/\hbar},&\mbox{if}\;0\leq
x\leq x_{\text{d}},\\\
0,&\mbox{if}\;x<0\;\text{and}\;x>x_{\text{d}},\end{cases}$ (37)

where $\bar{p}_{0x}$ and $E$ are the component $x$ of the momentum and the
energy of the system, respectively. Right after passing the slits, the $y$
component of the particle wave function reads

$\displaystyle\psi_{y}(y,0)=\frac{g_{-}\,\mbox{e}^{-i\bar{p}_{0y}(y+y_{\text{s}})/\hbar}+g_{+}\,\mbox{e}^{+i\bar{p}_{0y}(y-y_{\text{s}})/\hbar}}{\sqrt{2\left[1+\exp\left(-\tfrac{{y}_{\text{s}}^{2}}{2\sigma_{0}^{2}}+\frac{2\bar{p}_{0y}^{2}\sigma_{0}^{2}}{\hbar^{2}}\right)\right]}},$
(38)

where $g_{\pm}\equiv\sqrt{G_{y}(\pm{y}_{\text{s}},\sigma_{0})}$. In the $x$
direction, the Bohmian trajectory is no different from the classical one:
$x=x_{0}+\bar{p}_{0x}t/m$. Using the prescriptions (8) and (9), we obtain the
distributions

$\displaystyle{\wp_{x|t}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern
1.0mu}(x|t)=\begin{cases}x_{\text{d}}^{-1},&\mbox{if}\;0\leq x\leq
x_{\text{d}},\\\ 0,&\mbox{if}\;x<0\;\text{and}\;x>x_{\text{d}},\end{cases}$
(39)

and

$\displaystyle{\wp_{t|x}}^{\scriptscriptstyle\mkern-25.0mu\text{BM}\mkern
1.0mu}(t|x_{\text{d}})=\begin{cases}t_{\text{d}}^{-1},&\mbox{if}\;0\leq t\leq
t_{\text{d}},\\\ 0,&\mbox{if}\;t<0\;\text{and}\;t>t_{\text{d}}.\end{cases}$
(40)

where $t_{\text{d}}\equiv mx_{\text{d}}/|\bar{p}_{0x}|$ is the classical time
taken by a particle of mass $m$ that crosses the slit and moves straight with
momentum $\bar{p}_{0x}$ to reach the screen. Notice that, since the spatial
distribution is uniform in the interval $x\in[0,x_{\text{d}}]$, the flight-
time distribution is also uniform in the time interval
$t\in\left[0,t_{\text{d}}\right]$. Thus, by substituting the formulas above
into the prescription (2), we find

$\displaystyle{\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{BM}\mkern
6.0mu}(y|x_{\text{d}})=\frac{1}{t_{\text{d}}}\int_{0}^{t_{\text{d}}}\\!\\!dt\,|\psi_{y}(y,t)|^{2}.$
(41)

For added convenience, the integral will be carried out numerically.

In search of a counterpart to our result, we now supplement quantum mechanics
(QM) with a (semi)classical element recognized as suitable in far-field
regions. Specifically, we just substitute the classical time of flight,
$t_{\text{d}}=mx_{\text{d}}/|\bar{p}_{0x}|\equiv t_{\text{d}}(x_{\text{d}})$,
in the prescription (1), in order to obtain
${\wp_{y|x,t}}\big{(}y|x_{\text{d}},t_{\text{d}}(x_{\text{d}})\big{)}$, here
denoted as

$\displaystyle{\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{CT}\mkern
6.0mu}(y|x_{\text{d}})=|\psi_{y}\big{(}y,t_{\text{d}}(x_{\text{d}})\big{)}|^{2},$
(42)

with “CT” standing for classical time. Such a result is expected to be
appropriate solely for describing the far-field approximation scenario. As
done previously, we construct a similar probability distribution by counting
the number of Bohmian trajectories that cross the point $(x_{\text{d}},y)$
from an ensemble of $10^{4}$ trajectories with initial conditions
$(x_{0},p_{0})$ distributed according to
$x_{\text{d}}^{-1}|\psi_{y}(y_{0},0)|^{2}$. The result obtained from this
direct counting is depicted in Fig. 4 with blue dots. The results obtained
with our formalism, using BM as TET, are plotted as a solid black curve, while
the quantum-mechanical approach supplemented with a classical time of flight
is plotted as a solid red line.

Figure 4: Probability distributions of type $\wp_{y|x}$ as a function of $y$,
for the double-slit experiment, conditioned to a detection at the coordinate
$x=x_{\text{d}}$. The solid black curve represents the results from our
formalism ${\wp_{y|x}}^{\scriptscriptstyle\mkern-27.0mu\text{BM}\mkern
3.0mu}(y|x_{\text{d}})$, the red curve represents the results from quantum
mechanics supplemented with classical time
${\wp_{y|x}}^{\scriptscriptstyle\mkern-27.0mu\text{CT}\mkern
3.0mu}(y|x_{\text{d}})$, and the blue dots represent the results obtained from
the direct counting of an ensemble of $10^{4}$ Bohmian trajectories, with
initial conditions distributed according to
$x_{\text{d}}^{-1}|\psi_{y}(y_{0},0)|^{2}$. Panels (a) and (b) depict the
distributions for two sets of parameters:
${\\{\sigma_{0}=1,\bar{y}_{0}=0,\bar{p}_{0x}=1,\bar{p}_{0y}=2,m=0.5,x_{\text{d}}=1,y_{\text{s}}=2,\hbar=1\\}}$
and
${\\{\sigma_{0}=1,\bar{y}_{0}=0,\bar{p}_{0x}=1,\bar{p}_{0y}=2,m=0.5,x_{\text{d}}=3,y_{\text{s}}=2,\hbar=1\\}}$,
respectively. The only difference is the distance $x_{\text{d}}$ from the
detection screen to the slits. All physical quantities are in arbitrary units.

Although not always accentuated, the differences between the models
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{BM}\mkern
6.0mu}(y|x_{\text{d}})$ and
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{CT}\mkern
6.0mu}(y|x_{\text{d}})$ are perceptible and, in principle, measurable. In the
near-field regime [panel (a)], where the two Gaussian branches of the wave
function overlap significantly, both models predict the central interference
pattern, but with different visibilities. In the far-field regime [panel (b)],
the differences are remarkable and can be traced back to formulas (41) and
(42). When we take $x_{\text{d}}$ to be large, and thus the flight-time
$t_{\text{d}}=mx_{\text{d}}/|p_{0x}|$, the latter formula allows the branches
of the wave function to move far apart, reducing the overlap and preventing
the central interference pattern from manifesting itself. On the other hand,
formula (41) prescribes time averaging over several times, thus preserving the
central interference pattern, albeit with low visibility. In other words, from
an experimental standpoint,
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{CT}\mkern
6.0mu}(y|x_{\text{d}})$ corresponds to statistically post-selecting
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{BM}\mkern
6.0mu}(y|x_{\text{d}})$ for the specific time value $t_{\text{d}}$. As for the
counting result, it is noteworthy that it consistently corroborates our
statistical formalism, which is based on
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{BM}\mkern
6.0mu}(y|x_{\text{d}})$.

As a final remark, we note that the lower visibility exhibited by
${\wp_{y|x}}^{\scriptscriptstyle\mkern-30.0mu\text{BM}\mkern
6.0mu}(y|x_{\text{d}})$ in Fig. 4(a) might lead one to believe, when
confronted with experimental data of this form, that the system is subjected
to some source of decoherence. However, what is being shown here is that this
effect derives solely from flight-time fluctuations; no spurious,
uncontrollable degrees of freedom have been considered in our treatment. This
raises an intriguing question about low-visibility experimental patterns. If
one could guarantee the isolation of the system, then the low visibility could
be seen as corroborating our theory.

##  IV Final remarks

In this work, we demonstrate that trajectory-equipped theories have a
potential advantage over quantum mechanics in predicting the statistical
outcomes of experiments involving flight-time fluctuations. Specifically, we
show that by starting with the determinism hypothesis, a formalism can be
developed [see Eqs. (8), (9), (22), (24), and (26)] that yields predictions
for the flight-time probability distributions in emblematic experiments. From
a theoretical viewpoint, this result is significant because it raises the
perspective that Bohmian mechanics is not just an alternative interpretation
of the standard quantum formalism, but rather provides predictions unreachable
by the latter. More fundamentally, if confirmed empirically, these predictions
will provide support for the idea of an underlying trajectory-based
determinism spanning the quantum realm. On the other hand, if these
predictions are not confirmed, then Bohmian mechanics would be ruled out. We
hope that this work may stimulate researchers to test our results in the
laboratory, which, we believe, is feasible with current technology.

##  V ACKNOWLEDGMENTS

This study was financed in part by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. R.M.A. and
A.D.R. thank the financial support from the National Institute for Science and
Technology of Quantum Information (CNPq, INCT-IQ 465469/2014-0). R.M.A. also
thanks the Brazilian funding agency CNPq for Grant No. 305957/2023-6.

## References

  * [1] A. Zeilinger, R. Gähler, C. G. Shull, W. Treimer, and W. Mampe, Rev. Mod. Phys. 60, 1067 (1988).
  * [2] R. Gähler, and A. Zeilinger, Am. J. Phys. 59, 316-324 (1991).
  * [3] A. Tonomura, J. Endo, T. Matsuda, TT. Kawasaki, and H. Ezawa, Am. J. Phys. 57, 117-120 (1989).
  * [4] R. Bach, D. Pope, S.-H. Liou, and H. Batelaan, New J. Phys. 15, 033018 (2013).
  * [5] O. Carnal, and J. Mlynek, Phys. Rev. Lett. 66, 2689 (1991).
  * [6] C. Kurtsiefer, T. Pfau, and J. Mlynek, Nature 386, 150-153 (1997).
  * [7] C. Brand, S. Troyer, C. Knobloch, O. Cheshnovsky, and M. Arndt, Am. J. Phys. 89, 1132-1138 (2021).
  * [8] L. Maccone, and K. Sacha, Phys. Rev. Lett. 124, 110402 (2020).
  * [9] E. O. Dias, and F. Parisio, Phys. Rev. A 95, 032133 (2017).
  * [10] D. R. Terno, Phys. Rev. A 89, 042111 (2014).
  * [11] D. L. Sombillo, and E. A. Galapon, Ann. Phys. 364, 261-273 (2016).
  * [12] A. A. Rafsanjani, M. Kazemi, A. Bahrampour, and M. Golshani, Commun. Phys. 6, 195 (2023).
  * [13] N. Vona, G. Hinrichs, and D. Dürr, Phys. Rev. Lett. 111, 220404 (2013).
  * [14] S. Das, M. Nöth, and D. Dürr, Phys. Rev. A 99, 052124 (2019).
  * [15] D. Shucker, J. Funct. Anal. 38, 146-155 (1980).
  * [16] S. Wolf, and H. Helm, Phys. Rev. A 62, 043408 (2000).
  * [17] B. Korzh, et. al., Nat. Photonics 14, 250-255 (2020).
  * [18] S. Steinhauer, S. Gyger, and V. Zwiller, Appl. Phys. Lett 118, 100501 (2021).
  * [19] H. Azzouz, et. al., AIP Adv. 2, 032124 (2012).
  * [20] M. Rosticher, et. al., Appl. Phys. Lett. 97, 183106 (2010).
  * [21] F. Delgado, J. G. Muga, and G. García-Calderón, Phys. Rev. A 74, 062102 (2006).
  * [22] D. Bohm, Phys. Rev. 85, 166 (1952).
  * [23] D. Bohm, Phys. Rev. 85, 180 (1952).
  * [24] C. R. Leavens, Phys. Rev. A 58, 840 (1998).
  * [25] J. G. Muga, and C. R. Leavens, Phys. Rep. 338, 353-438 (2000).
  * [26] S. Das, and D. Dürr, Scientific Rep. 9, 2242 (2019).
  * [27] E. Nelson, Phys. Rev. 150, 1079 (1966).
  * [28] M. J. W. Hall, D.-A. Deckert, and H. M. Wiseman, Phys. Rev. X 4, 041013 (2014).
  * [29] F. McLafferty, The J. Chem. Phys. 83, 5043-5045 (1985).
  * [30] H. Nitta, and T. Kudo, Phys. Rev. A 77, 014102 (2008).
  * [31] M. J. Kazemi, and V. Hosseinzadeh, Phys. Rev. A 107, 012223 (2023).
  * [32] A. Ayatollah Rafsanjani, M. J. Kazemi, V. Hosseinzadeh, and M. Golshani, Scientific Rep. 14, 3615, (2024).
  * [33] J. Kijowski, Rep. Math. Phys. 6, 361-386 (1974).
  * [34] I. L. Egusquiza, and J. G. Muga, Phys. Rev. A 61, 012104 (1999).
  * [35] N. Page, and K. Wootters, Phys. Rev. D 27, 2885 (1983).
  * [36] J. S. Briggs, and J. M. Rost, Eur. Phys. J. D 10, 311-318 (2000).
  * [37] Eduardo O. Dias, Phys. Rev. A 103, 012219 (2021).
  * [38] G. Sebastian, and J. M. Rost, Phys. Rev. Lett. 131, 140202 (2023).
  * [39] R. Ximenes, F. Parisio, and E. O. Dias, Phys. Rev. A 98, 032105 (2018).
  * [40] R. Werner, Ann. l’I.H.P. Phys. Théorique 47, 429-440 (1987).
  * [41] R. Tumulka, Ann. Phys. 442, 168910 (2022).
  * [42] R. Tumulka, Phys. Rev. A 106, 042220 (2022).
  * [43] V. Dubey, C. Bernardin, and A. Dhar, Phys. Rev. A 103, 032221 (2021).
  * [44] A. Marchewka, and Z. Schuss, Phys. Rev. A 63, 032108 (2001).
  * [45] A. Marchewka, and Z. Schuss, Phys. Lett. A 240, 177-184 (1998).
  * [46] A. Marchewka, and Z. Schuss, Phys. Rev. A 61, 052107 (2000).
  * [47] T. Jurić, and H. Nikolić, Eur. Phys. J. Plus 137, 631 (2022).
  * [48] S. Roncallo, K. Sacha, and L. Maccone, Quantum 7, 968 (2023).
  * [49] J. A. Damborenea, I. L. Egusquiza, G. C. Hegerfeldt and J. C. Muga, Phys. Rev. A 66, 052104 (2022).
  * [50] J. C. Muga, S. Brouard, and D. Macías, Ann. Phys. 240, 351-366 (1995).
  * [51] J. J. Halliwell, Phys. Rev. A 240, 237-242 (1995).
  * [52] A. J. Bracken, and G. F. Melloy, J. Phys. A: Math. Gen. 27, 2197 (1994).
  * [53] H. Hofmann, Phys. Rev. A 96, 020101 (2017).
  * [54] K. Husimi, Proc. Phys. Math. Soc. Jpn. 22, 264-314 (1940).
  * [55] W. Greiner, Classical mechanics: Systems of particles and Hamiltonian dynamics (Springer-Verlag, Berlin, 2010).


[a]A. Tewsley-Booth

# Magnetic Field Analysis for Fermilab Muon $g-2$

_for the Muon $g-2$ collaboration_

###### Abstract

The Muon g-2 experiment E989 at Fermilab measures the anomalous magnetic
moment of the muon $a_{\mu}$ with improved precision compared to the
Brookhaven (E821) experiments. The Brookhaven results are in tension with the
Standard Model by more than $3\sigma$. The determination of $a_{\mu}$ requires
the measurement of both the muon anomaly frequency, $\omega_{a}$, and the
magnetic field, $\mathbf{B}$, that confines muons in a storage ring. The field
is monitored by a set of coordinated nuclear magnetic resonance (NMR)
measurements. NMR probes at fixed locations above and below the storage region
constantly monitor the field. An in-vacuum trolley equipped with 17 NMR probes
maps the muon storage region, and a special water-based NMR probe provides the
calibration for the trolley probes. This presentation focuses on the
determination of the time-dependent field maps from combining the fixed probe
measurements and the trolley maps. The field maps are combined with the muon
distribution to derive the average field observed by the muons during the
measurement. These proceedings will cover the analysis from the first data
run.

##  1 Introduction

Our collaboration reported a new measurement of the muon magnetic anomaly,
$a_{\mu}=\frac{1}{2}(g_{\mu}-2)$, based on our Run-1 data set, which was
collected between March and July of 2018. The result is

$a_{\mu}(\mathrm{FNAL})=116592040(54)\times 10^{-11},$ (1)

measured to a precision of 460 ppb, and is in good agreement with the previous
measurement from Brookhaven [2]. A set of four companion papers cover the
final result [3], the anomalous spin precession frequency [4], the magnetic
field measurements [5], and the beam dynamics corrections [6]. These
proceedings will broadly cover the magnetic field measurement and analysis for
the Run-1 data set.111A similar talk was given by the speaker at the 2021
DPyC-SMF conference.

The spin of particle with a magnetic moment precesses about an external non-
parallel magnetic field at the Larmor frequency [7],

$\omega=\frac{gq}{2m}B$ (2)

. If the particle’s charge, $q$, and its mass, $m$, are known, simultaneous
measurement of the spin precession frequency $\omega$ and the magnetic field
$B$ would allow the determination of the g-factor of the particle. The muon’s
magnetic moment is slightly greater than two due to interactions with the
vacuum. This shift, the magnetic anomaly $a=\frac{g-2}{2}$, is effected by
interactions with virtual particles in the vacuum. Therefore, a measurement of
the magnetic anomaly can be used to probe properties of these interactions.
Current $g-2$ theory takes into account QED, electroweak, and QCD. The QED
terms dominate the value, but the QCD terms dominate the uncertainty [8].

In a constant magnetic field, the relative angle between the momentum and spin
vectors of the muon evolves over time. The rate of change of this angle,
called the anomalous precession frequency $\omega_{a}$, is (nearly)
proportional to the magnetic anomaly, $a_{\mu}$.

$\vec{\omega}_{a}=\vec{\omega}_{s}-\vec{\omega}_{c}\approx-\frac{q}{m}a_{\mu}\vec{B},$
(3)

where $\omega_{s}$ is the spin precession frequency and $\omega_{c}$ is the
cyclotron (momentum precession) frequency. Simultaneous measurement of
$\omega_{a}$ and the magnetic field enable determination of $a_{\mu}$. In this
experiment, the field is measured as $\tilde{\omega}^{\prime}_{p}$, the Larmor
frequency of the proton in water in the same magnetic field, weighted by the
muon distribution.

The muon-averaged field, $\tilde{\omega}^{\prime}_{p}$, can be laid out
schematically as

$\tilde{\omega}^{\prime}_{p}=f_{\text{field}}~{}\left\langle\omega_{\mathrm{p,meas}}\bigotimes\rho_{\mu}\right\rangle~{}(1+B_{q}+B_{k}).$
(4)

Here, $f_{\mathrm{field}}$ is the absolute calibration of the field that takes
$\omega_{\text{p, meas}}$ to $\omega^{\prime}_{p}$, the equivalent Larmor
frequency of the proton in water. The term
$\left\langle\omega_{\mathrm{p,meas}}\bigotimes\rho_{\mu}\right\rangle$ is the
average measured field weighted by the muon distribution $\rho_{\mu}$. The
terms $B_{q}$ and $B_{k}$ are the effects of fast transient magnetic fields
that cannot be measured by the field systems that measure $\omega_{\text{p,
meas}}$.

##  2 Measuring $\omega_{p}$

The magnetic field experienced by the muons is measured using several
different magnetometer systems. The most important of these are the
calibration probe, trolley, and fixed probe systems. The primary field mapping
system, the trolley, is an array of magnetometers located in the vacuum
chamber. The trolley is pulled by cable around the muon storage region,
measuring the field as it goes. Each of these systems uses NMR magnetometry,
and together they form a calibration chain that allows us to measure the
absolute field to a precision of 114 ppb (56 ppb from the calibration,
measurements, analysis, and averaging; and 99 ppb from the effects of fast
transient fields).

Figure 1: The trolley generates precise field maps as a function of azimuthal
position as it travels around the ring. Between trolley runs, these maps are
interpolated using measurements from the fixed probes that track the evolution
of the field.

The calibration chain begins at Argonne National Laboratory in a precision MRI
magnet. There, the calibration probe is cross-checked with a He-3 NMR probe
[9]. The two probes were found to be in excellent agreement. The calibration
probe is then transferred to the storage ring at FNAL, where the calibration
from the calibration probe is transferred to each of the trolley probes. Then,
as the trolley moves about the storage ring, it can transfer the calibration
to the array of fixed probes.

Figure 2: The relative positions of the fixed and trolley probes in an
azimuthal slice of the ring. In this coordinate system, the muon’s magic
radius is at $(0,0)$.

After calibration, the primary systems used to map the magnetic field in the
storage ring are the trolley and the fixed probes. The trolley makes its
measurements in the same volume that the muons fill. It has 17 NMR probes and
takes measurements in about 4000 azimuthal locations around the ring, as shown
in Figure 1. However, the trolley cannot be operated during muon fills because
it blocks their path, so it is pulled out of the way into its garage most of
the time, only mapping the field about every three days. On the contrary, the
fixed probes only measure at 72 azimuthal locations and are physically located
outside the vacuum chambers, but they can continue measuring the field during
muon fills, providing information about how the field evolves between trolley
runs. Figure 2 shows the relative locations of the trolley and fixed probes in
an azimuthal slice of the ring. The two sources are combined in the analysis,
with the fixed probe data being used to interpolate the field map between the
trolley runs. For Run-1, this procedure, from measurement through analysis,
accounts for 56 ppb of the total uncertainty.

The NMR magnetometers are good for measuring quasi-static fields, but
transients with characteristic times faster than about a second require
special systematic studies to measure. There are two primary sources of
transient magnetic fields in the storage ring: the electrostatic quadrupoles
and the faster kickers. The ESQ plates are charged to high voltage at each
muon injection to provide vertical focusing on the beam. Charging the plates
induces vibrations, which generate a magnetic field that perturbs the muons.
For the Run-1 analysis, the effect of this transient field was found to be -17
ppb with an uncertainty of 92 ppb, making it the dominant systematic effect on
the final result. Since the Run-1 publication, we have completed a more in-
depth systematic study and expect this uncertainty to be reduced in future
results. Furthermore, beginning in Run-5, we have changed the ESQ charging
procedure to reduce the induced vibrations that cause the transient.

Muon OrbitStraw TrackersCalorimeterHigh Momentum Decay PositronLow Momentum
Decay Positron Figure 3: Decay positrons travel through the straw trackers.
Their paths can be traced back using knowledge of the magnetic field to their
decay vertices, allowing for a calculation of the muon distribution.

The source of the other transient is the kicker system, which uses a fast
magnetic field at the beginning of each muon injection to kick the muons onto
their ideal orbit. The kick induces eddy currents that perturb the field as
they decay. This effect is measured using a Faraday magnetometer that can
measure the field at the nanosecond level. The average effect on the muons was
found to be -27 ppb, with an uncertainty of 37 ppb. This measurement has also
been refined and repeated since the Run-1 publication to reduce its associated
uncertainty. Table 1 details the corrections and uncertainties for the quasi-
static and transient dynamics of the field, weighted by the muon distribution.

In order to calculate the average magnetic field experienced by the muons, we
need to know both the magnetic field and the muon distribution in the storage
ring as a function of position and time. The muon distribution is measured by
the straw trackers at two azimuthal locations. Those two distributions are
then used to extrapolate the distribution around the whole ring by combining
the measurement with beam dynamics simulations.

The straw trackers are formed of layers of overlapping straws filled with gas
that is ionized as positrons travel through the device towards the
calorimeters. By measuring the positions where the ionization occurs in the
straws, we can find the tracks the positrons took from their decay positions
(see Figure 3). Extrapolating backwards through the magnetic field, these
tracks are used to determine the decay vertices of the positrons, which are
used as a proxy for the muon distribution at the two azimuthal locations of
the straw trackers. We use the distributions averaged over several hours of
data collection as a weighting function when we average the magnetic field’s
non-uniformity.

Quantity | Correction Terms | Uncertainty
---|---|---
| (ppb) | (ppb)
$\tilde{\omega}^{\prime}_{p}$ | – | 56
$B_{k}$ | -27 | 37
$B_{q}$ | -17 | 92
Totals | -44 | 114
Table 1: The corrections and uncertainties to the field terms in Equation 4.

##  3 Acknowledgments

We thank the Fermilab management and staff for their strong support of this
experiment, as well as the tremendous support from our university and national
laboratory engineers, technicians, and workshops. The Muon $g-2$ Experiment
was performed at the Fermi National Accelerator Laboratory, a U.S. Department
of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi
Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359.
Additional support for the experiment was provided by the Department of Energy
offices of HEP and NP (USA), the National Science Foundation (USA), the
Istituto Nazionale di Fisica Nucleare (Italy), the Science and Technology
Facilities Council (UK), the Royal Society (UK), the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie grant
agreements No. 690835, No. 734303, the National Natural Science Foundation of
China (Grant No. 11975153, 12075151), MSIP, NRF and IBS-R017-D1 (Republic of
Korea), the German Research Foundation (DFG) through the Cluster of Excellence
PRISMA+ (EXC 2118/1, Project ID 39083149).

## References

  * [1]
  * [2] G. Bennett _et. al._ , Phys. Rev. D 73, 7 (2006).
  * [3] B. Abi _et. al._ , Phys. Rev. Lett. 126, 14 (2021).
  * [4] T. Albahri _et. al._ , Phys. Rev. D 103, 7 (2021).
  * [5] T. Albahri _et. al._ , Phys. Rev. A 103, 4 (2021).
  * [6] T. Albahri _et. al._ , Phys. Rev. Accel. Beams 24, 4 (2021).
  * [7] J. D. Jackson. Classical Electrodynamics, Third Edition. Wiley (1998).
  * [8] T. Aoyama _et. al._ , Phys. Reports, 887, (2020).
  * [9] M. Farooq _et. al._ , Phys. Rev. Lett. 124, 22, (2020).


Ginwidth=,height=,keepaspectratio

# scida: scalable analysis for scientific big data

Chris Byrohl  Heidelberg University, Institute for Theoretical Astronomy,
Albert-Ueberle-Str. 2, 69120 Heideberg, Germany  0000-0002-0885-8090 Dylan
Nelson  Heidelberg University, Institute for Theoretical Astronomy, Albert-
Ueberle-Str. 2, 69120 Heideberg, Germany  0000-0001-8421-5890

(15 Sep 2023)

## Summary

scida (https://scida.io) is a Python package for reading and analyzing large
scientific data sets. Data access is provided through a hierarchical
dictionary-like data structure after a simple load() function. Using the dask
library for scalable, parallel and out-of-core computation (Dask Development
Team, 2016), all computation requests from a user session are first collected
in a task graph. Arbitrary custom analysis, as well as all available dask
(array) operations, can be performed. The subsequent computation is executed
only upon request, on a target resource (e.g. a HPC cluster, see Figure 1).

Figure 1: Schematic of the workflow. In a user session, a recipe (i.e.
sequence of analysis operations) for desired data product can be built by
consecutive chaining of operations, which are internally represented by dask
task graphs. Calculation is triggered by the compute() command, evaluating the
graph on a target resource. The result, much smaller than the original data,
is sent back to the user session for further analysis/plotting.

## Features

Scida begins by providing a clean, abstract, dictionary-like interface to the
underlying data, regardless of its file format or structure on disk. Physical
units are automatically inferred, when possible, and attached to the dataset.
Symbolic and automatic unit conversion is provided by the pint package (Grecco
& others, 2012). This metadata can also be specified by the user via
customizable configuration files.

Scida attempts to automatically determine the type of dataset loaded. At the
file level, it currently supports zarr, single/multi-file HDF5 and FITS files.
The analysis functionality then available depends on the dataset type. For
example:

  * •

Datasets which pair scattered pointsets with cataloged results of clustering
algorithms have a natural association between such groups and their member
points. Cosmological simulations in astrophysics provide one such motivating
example, where halos and galaxies are identified in particle data. Scida then
supports broadcasting of analysis and reduction operations over all groups for
the associated particles.

  * •

Datasets containing spatial coordinates in one, two, or three dimensions can
be queried for spatial subsets.

  * •

Datasets which are comprised of series, e.g. snapshots at different times, or
parameter variation suites, are automatically inferred and supported.

The dataset-dependent functionality of scida is handled with a flexible and
extensible mixin architecture.

## Statement of need

Today, scientific “big data” is petabyte-scale. This makes traditional
analysis by a researcher on their personal computer prohibitive. Writing
analysis code which is distributed or out-of-core is complex, and is not the
main focus of scientists. The need for significant data management creates (i)
a barrier for new researchers, (ii) a substantial time commitment apart from
the science itself, and (iii) an increased risk of errors and difficulties in
reproducibility due to higher code complexity, while (iv) workflows are often
not easily transferable to other datasets, nor (v) scalable and transferable
to changing computing resources.

scida solves these problems by providing a simple interface to large datasets,
hiding the complexity of the underlying data format and file structures, and
transparently handling the parallelization of analysis computations. This is
facilitated by the dask library, which naturally separates the definition of a
given computation from its execution.

Initial support in scida is focused on astrophysical simulations and
observations, but the package is designed to be easily extensible to other
scientific domains. Existing analysis frameworks for astrophysical simulations
include python packages such as yt (Turk et al., 2011), pynbody (Pontzen et
al., 2013), pygad (Röttgers, 2018), nbodykit (Hand et al., 2018) and
swiftsimio (Borrow & Kelly, 2021). None utilize the graph-based distributed
analysis framework of dask. Often, existing analysis packages rely on the
explicit loading of entire datasets into main memory. However, this approach
is not transparently scalable to large data sets, and requires the user to
explicitly manage the data chunks in custom analysis routines.

By providing a flexible interface to the dask library to handle large data
sets in a scalable fashion, users can also leverage dask functionality and
dask-based libraries such as dask-image (Kirkham & others, 2018) and
datashader (Bednar et al., 2022).

scida was first utilized in Byrohl & Nelson (2023) for the analysis of
cosmological and radiative transfer simulations, particularly the reduced chi-
squared analysis exploring thousands of terabyte-sized models.

## Target Audience

scida aims to simplify access to large scientific data sets. It lowers the
barrier of entry for researchers to ask complex questions of big data. As
such, the scida package is targeted at researchers with large data analysis
problems. Its clean interface is appropriate for users with no prior
experience with big data or distributed data analysis, as well as those who
specifically want to leverage dask to make their workflows easier to read and
scalable.

Domain specific analysis routines can be implemented on top of scida. Initial
“out of the box” data support is currently focused on astrophysical data sets,
but scida aims to support other scientific domains as well, where similar
solutions will be beneficial.

## Acknowledgements

CB and DN acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG)
through an Emmy Noether Research Group (grant number NE 2441/1-1).

## References

reBednar, J. A., Crail, J., Crist-Harif, J., Rudiger, P., Brener, G., B, C.,
Thomas, I., Mease, J., Signell, J., Liquet, M., Stevens, J.-L., Collins, B.,
Thorve, A., Bird, S., thuydotm, esc, kbowen, Abdennur, N., Smirnov, O., …
Bourbeau, J. (2022). _Holoviz/datashader: Version 0.14.3_ (Version v0.14.3).
Zenodo. https://doi.org/10.5281/zenodo.7331952

preBorrow, J., & Kelly, A. J. (2021). _Projecting SPH particles in adaptive
environments_. https://arxiv.org/abs/2106.05281

preByrohl, C., & Nelson, D. (2023). The cosmic web in Lyman-alpha emission.
_Monthly Notices of the Royal Astronomical Society_ , _523_ , 5248–5273.
https://doi.org/10.1093/mnras/stad1779

preDask Development Team. (2016). _Dask: Library for dynamic task scheduling_.

preGrecco, H., & others. (2012). Pint: Makes units easy. In _GitHub
repository_. GitHub. https://github.com/hgrecco/pint

preHand, N., Feng, Y., Beutler, F., Li, Y., Modi, C., Seljak, U., & Slepian,
Z. (2018). nbodykit: An Open-source, Massively Parallel Toolkit for Large-
scale Structure. _156_(4), 160. https://doi.org/10.3847/1538-3881/aadae0

preKirkham, J., & others. (2018). Dask-image: Distributed image processing.
_GitHub Repository_. https://github.com/dask/dask-image

prePontzen, A., Roškar, R., Stinson, G. S., Woods, R., Reed, D. M., Coles, J.,
& Quinn, T. R. (2013). _pynbody: Astrophysics Simulation Analysis for Python_.

preRöttgers, B. (2018). _pygad: Analyzing Gadget Simulations with Python_ (p.
ascl:1811.014). Astrophysics Source Code Library, record ascl:1811.014.

preTurk, M. J., Smith, B. D., Oishi, J. S., Skory, S., Skillman, S. W., Abel,
T., & Norman, M. L. (2011). yt: A Multi-code Analysis Toolkit for
Astrophysical Simulation Data. _The Astrophysical Journal Supplement Series_ ,
_192_ , 9. https://doi.org/10.1088/0067-0049/192/1/9

p