Deterministic and randomized motions in single-well potentials

TL;DR

This paper compares long-term behaviors of deterministic and randomized motions in single-well potentials, revealing that in the long run, probability distributions become insensitive to velocity reversals under certain conditions.

Contribution

It introduces a specific energy-conserving randomization method involving velocity reversals and analyzes its effects on long-term probability distributions in 1D and 2D systems.

Findings

Distributions become insensitive to velocity reversals in 1D long-term limit.

Additional randomization of initial conditions affects distributions.

In 2D, distributions are also insensitive to velocity reversals asymptotically.

Abstract

Newtonian, undamped motion in single-well potentials belong to a class of well-studied conservative systems. Here, we investigate and compare long-time properties of fully deterministic motions in single-well potentials with analogous randomized systems. We consider a special type of energy-conserving randomization process: the deterministic motion is interrupted by hard velocity reversals $\vec{v}(t_i)\to-\vec{v}(t_i)$ at random time instants $t_i$. In the 1D case, for fixed initial conditions, the differences in probability distributions disappear in the long-time limit making asymptotic densities insensitive to the selection of random time instants when velocity is reversed. Substantially different probability distributions can be obtained, for instance, through the additional randomization of initial conditions. Analogously, in 2D setups, the probability distributions asymptotically are insensitive to velocity reversals.

Equations

$$m\frac{d^{2}x(t)}{d{t}^2}=-V^{\prime }(x)$$

$$V(x)=\kappa |x{|}^n(\kappa >0,n>0)$$

$$T=\frac{2}{n}\sqrt{\frac{2\pi m}{E}}{\left[\frac{E}{\kappa }\right]}^{\frac{1}{n}}\frac{\Gamma \left(\frac{1}{n}\right)}{\Gamma \left(\frac{1}{2}+\frac{1}{n}\right)},$$

$$E=\frac{1}{2}m{v}^2+V(x)=\frac{1}{2}m{v}^2+\kappa |x{|}^n.$$

$$p(x)=\frac{2}{T}{\left[\frac{2}{m}\left(E-\kappa |x{|}^n\right)\right]}^{-\frac{1}{2}}.$$

$$p(v)=\frac{2m}{T\kappa n}{\left[\frac{1}{\kappa }\left(E-\frac{1}{2}m{v}^2\right)\right]}^{\frac{1}{n}-1}.$$

$$p(x)=\int p(x|\lambda )p(\lambda )d\lambda$$

$$p(v)=\int p(v|\lambda )p(\lambda )d\lambda ,$$

$$m\frac{d^2\vec{r}(t)}{d{t}^2}=-\nabla V(r).$$

$$p_x(x)=\frac{1}{\pi \sqrt{A_x^2-x^2}},$$

$$p_v(v_x)=\frac{1}{\pi \sqrt{(A_x\omega {)}^2-v_x^2}}.$$

$$p(r)\propto \frac{\sqrt{m}}{2\left[\sqrt{E-U_{\mathrm{eff}}(r)}\right]},$$

$$U_{\mathrm{eff}}(r)=\frac{L^2}{2m{r}^2}+\kappa r^n,$$

$$\varphi (k)=\exp [-|k{|}^{\alpha }].$$

$$p(x)dx$$

$$T=\frac{2}{n}\sqrt{\frac{2\pi m}{E}}{\left[\frac{E}{\kappa }\right]}^{\frac{1}{n}}\frac{\Gamma \left(\frac{1}{n}\right)}{\Gamma \left(\frac{1}{2}+\frac{1}{n}\right)}.$$

$$E=\frac{1}{2}m{v}^2+V(x)=\frac{1}{2}m{v}^2+\kappa |x{|}^n$$

$$v=\sqrt{2\left[E-\kappa |x{|}^n\right]}.$$

$$p(x)=\frac{2}{T}\frac{1}{\sqrt{2\left[E-\kappa |x{|}^n\right]}}.$$

$$p(v)=p(x(v))\left|\frac{dx}{dv}\right|.$$

$$x={\left[\frac{1}{\kappa }\left(E-\frac{1}{2}m{v}^2\right)\right]}^{\frac{1}{n}}$$

$$\frac{dx}{dv}=-\frac{m}{\kappa n}{\left[\frac{1}{\kappa }\left(E-\frac{1}{2}m{v}^2\right)\right]}^{\frac{1}{n}-1}v.$$

$$p(v)=\frac{2m}{T\kappa n}{\left[\frac{1}{\kappa }\left(E-\frac{1}{2}m{v}^2\right)\right]}^{\frac{1}{n}-1}.$$

$$\tilde{p}(E,x)=\frac{En\sqrt{\frac{m}{E}}\Gamma \left(\frac{1}{2}+\frac{1}{n}\right){\left(\frac{E}{\kappa }\right)}^{-1/n}{\left(\frac{E{|x|}^{-n}}{\kappa }\right)}^{\frac{1}{n}-\frac{1}{2}}{}_2{F}_1\left(\frac{1}{2},\frac{1}{n}-\frac{1}{2};\frac{3}{2};1-\frac{E{|x|}^{-n}}{\kappa }\right)\sqrt{\frac{E-\kappa {|x|}^n}{m}}}{\sqrt{\pi }\Gamma \left(\frac{1}{n}\right)}$$

$$\tilde{p}(E,v)=$$

$$p(x)=\frac{\tilde{p}(E_2,x)-\tilde{p}(E_1,x)}{E_2-E_1}$$

$$p(v)=\frac{\tilde{p}(E_2,v)-\tilde{p}(E_1,v)}{E_2-E_1}.$$

$$\left\{\begin{array}[]{lcl}m\ddot{x}(t)&=&-2\kappa x(t)\\ m\ddot{y}(t)&=&-2\kappa y(t)\\ \end{array}\right.,$$

$$x(t)=A_x\cos (\omega t+{\Delta }_x),$$

$$p_x(x)dx=\frac{2}{T}\frac{dx}{\dot{x}}.$$

Full text

1 Introduction and motivation

Newtonian dynamics provides an accurate, deterministic description of macroscopic systems. The prerequisite for its application is the knowledge of the initial conditions. In the absence of dissipation, the Newtonian dynamics in static potentials is conservative, i.e., mechanical energy is constant. Consequently, the motion in the phase space continues along a constant energy curve. Variable, stochastic environment can be introduced by time-dependent parameters or by the introduction of noise [1, 2]. Noise is usually considered to be an approximate description of interactions between the studied particle/body and its environment. The action of additive noise, which is not counterbalanced by some other mechanism, e.g., dissipation, breaks both mean and instantaneous energy conservation [3, 4]. Nevertheless, a special class of disturbances in dynamical systems can still preserve the energy of the system, see Dybiec et al. [5]. For instance, the disturbance which reverses the velocity vector at random instants of time does not break the energy conservation at all. Presence of the stochastic component in the system dynamics transforms state variables into random variables. Properties of perturbed systems can be studied by means of statistical physics. Analogously, fully deterministic systems can also be studied via notions of probability theory. The long observation of deterministic trajectories allows one to estimate the probability density of finding a system in a given state [6].

Random perturbations in dynamical systems play a crucial role in noise-driven systems. Among noise-induced effects there are well-known and celebrated effects of stochastic resonance [2, 7], resonant activation [8, 9, 10], stochastic synchronization [11] and direct transport (ratcheting effect) [12]. The presence of noise can also induce dynamical multimodality [13, 14, 15] and dynamical hysteresis [16, 17]. All these effects originate in or owe its efficiency due to the combined action of deterministic and random forces. For instance, in the stochastic resonance, the combined action of the deterministic perturbation and fine-tuned noise amplify the system response by making the random process of passages over the potential barrier optimally synchronized with the periodic modulation of the potential barrier.

Current investigations are motivated by the Lévy walk model [18]. The Lévy walk is a dynamical process with simple spatiotemporal coupling. The random walker moves with a constant velocity $v$ for a random independent and identically distributed times $\tau_{i}$ ($\tau_{i}>0$) followed by a velocity flip, i.e., $v\to-v$. The time $\tau$ determines the flight time and the distance traveled in between velocity reversals. The mechanical energy of the Lévy walker is constant and equal to kinetic energy. Lévy walks were successfully applied to a plenitude of models including, among others, dynamics of tracer particles in weakly chaotic systems [19, 20, 21], cold atom systems [22], random search strategies [23, 24] and intracellular motion [25, 26].

In Dybiec et al. [5] the extension of the Lévy walk scenario including deterministic forces has been suggested. For the deterministic motion in a single-well potential velocity continuously changes due to the deterministic $-V^{\prime}(x)$ force. Therefore, in 1D, the absolute value of velocity $|v(t)|$ decreases as the particle moves towards large $|x|$ and at large $|x|$ the motion is smoothly reversed when $|v(t)|=0$. In addition to deterministic soft velocity reversals, hard velocity reversals at random time instants $t_{i}$ have been introduced. Hard velocity reversals are responsible for velocity flips $\vec{v}(t_{i})\to-\vec{v}(t_{i})$.

Here, we summarize some of the findings and further extend the model of Ref. [5]. Using the specific type of a multiplicative disturbance which reverses the velocity vector at random instants of time the deterministic system is randomized. The process of velocity reversals can be also interpreted as a dichotomous (assuming two values $\pm 1$), multiplicative perturbation $\eta(t)$ to the velocity $\vec{v}(t)\to\eta(t)\vec{v}(t)$. On the one hand, this process does not break the energy conservation. On the other hand, it introduces a discontinuity in the phase space because of discontinuous velocity changes. Nevertheless, the system stays on the same orbit, and asymptotically its statistical properties are the same as properties derived from the very long observation of the fully deterministic motion. Significant changes in its statistical properties are observable due to a different type of randomizations, e.g., randomization in its initial conditions. Randomization of initial conditions makes the total energy random variable and produce densities which significantly deviate from $u$-shaped densities appearing for fixed initial conditions accompanied with hard velocity reversals [5]. Within the studied model, the main source of the randomization is due to hard velocity reversals. Therefore, the model resembles the model of hard spheres [27] in the square-well potential which in 2D, two-particle, molecular dynamics ensemble is able to produce similar probability density functions for the velocity like the model studied within the current manuscript. Nevertheless, both models bear fundamental differences: the setup studied here is a single particle and non-zero potential model of randomized motion while the one considered in [27] is a multi-particle 2D and 3D model of the hard sphere gas.

The model under consideration is presented in the next section (Sec. 2). Section Results (Sec. 3) presents the outcome of our analysis. Starting from the conclusions from Dybiec et. al. [5] about hard velocity reversals we provide universal, analytical derivation of general formulas which, for simplest cases of single-well potentials, were constructed semi-analytically and numerically in [5]. Next, we ask an ancillary question about the long-time behavior of a class of fully deterministic systems in which initial conditions are randomized (Section 3.1). Finally, the higher dimensional case (2D) of quasi-periodic dynamics in 2D single-well potentials is analyzed with and without hard velocity reversals (Section 3.2). The manuscript is closed with Summary and Conclusions (Sec. 4). More technical derivations are moved into the appendices (A – D).

2 Model

The Newton’s equation

[TABLE]

with

[TABLE]

describes a periodic motion [28] with the period $T$ given by

[TABLE]

where $\Gamma(\dots)$ is the Euler Gamma function. Due to the absence of damping the system described by Eq. (1) is conservative. Its total energy $E$ is equal to

[TABLE]

From a single trajectory, after a long observation time, it is possible to calculate the probability of observing the particle in the vicinity of the point $x$

[TABLE]

Analogously, the probability of recording velocity $v=\dot{x}$ reads

[TABLE]

In Eqs. (5) and (6) $T$ stands for the period of the motion which is given by Eq. (3). The total energy $E$ is fixed by the initial condition. Densities given by Eqs. (5) and (6) are similar to the arcsine distribution with modes located at $V(x)=E$ for $p(x)$ or $mv^{2}/2=E$ for $p(v)$, see Fig. 1. For selected values of $n$ Eqs. (5) and (6) were derived in [5] using semi-analytical and numerical methods. The detailed, universal, analytical derivation of general formulas (5) and (6) is presented in A.

In Dybiec et al. [5], the 1D deterministic motion governed by Eq. (1) with fixed initial condition is perturbed by hard velocity reversals which do not break the energy conservation. Such hard velocity reversals at random time instants are separated by a random time $\tau$. More precisely, the sequence of waiting times $\tau_{i}$ follows a one-sided probability density $p(\tau)$. The particle moves according to Eq. (1) as long as $t\leqslant\tau_{1}$. Then, the velocity is reversed $v(\tau_{1})\to-v(\tau_{1})$ and the particle moves deterministically till the next velocity reversal, i.e., while $\tau_{1}<t\leqslant\tau_{1}+\tau_{2}$. The next velocity reversal is at time $\tau_{1}+\tau_{2}$, i.e., $v(\tau_{1}+\tau_{2})\to-v(\tau_{1}+\tau_{2})$. The whole procedure is repeated iteratively. At time instants $\tau_{1}$, $\tau_{1}+\tau_{2},\dots$ only the velocity is reversed while the position is unaffected. Therefore, we solve Eq. (1) in time intervals $0\leqslant t\leqslant\tau_{1},\tau_{1}\leqslant t\leqslant\tau_{1}+\tau_{2},\dots$ with the initial conditions determined by the final position and the reversed final velocity from a preceding interval.

3 Results

First, we reexamine the periodic motion in 1D single-well potentials (Sec. 3.1). Next, we switch to the 2D single-well potentials (Sec. 3.2). We study similarities and differences between 1D and 2D models with particular attention to the consequences of the fact that orbits in 2D single-well potentials do not need to be closed.

3.1 1D randomized motions

As it was shown in Dybiec et al. [5] asymptotic $p(v)$ and $p(x)$ densities are robust to the exact shape of $p(\tau)$. At the same time, transient behavior can display some sensitivity to the $p(\tau)$. At short times randomized motion following Eq. (1) produces stationary states which can be decorated by fingerprints of initial conditions. This effect is the most persistent for heavy-tailed distributions of reversal times which are characterized by the diverging mean. Regardless of the waiting time distribution $p(\tau)$ peaks corresponding to the initial condition decay over time.

Asymptotically, hard velocity reversals at random time instants produce the same position and velocity distributions regardless of $p(\tau)$, i.e., asymptotics densities are given by Eqs. (5) and (6) with appropriately adjusted parameters. Consequently, we now have three scenarios which end up with the same $p(x)$ and $p(v)$ distributions, namely when we build those distributions out of:

1. (i)

very long observation of a single deterministic trajectory with a fixed initial condition, 2. (ii)

ensemble averaging of final velocities and positions with random initial conditions uniformly distributed over the fixed energy orbit, 3. (iii)

ensemble averaging of final velocities and positions with fixed initial condition accompanied with random velocity reversals.

The scenario (i) corresponds to the time averaging over a single very long trajectory, which evolves according to the Newton equation of motion. For the very long observation time, with the fixed frequency, we sample $x$ and $v$ which are used to estimate $p(x)$ and $p(v)$ densities. The protocol (ii) is the ensemble averaging of trajectories with the fixed total energy $E$ but different initial velocities and positions. More precisely, in (ii) we generate a large number of starting points (initial conditions) which are uniform on the fixed energy curve. For each initial condition we evolve $x(t)$ and $v(t)$ according to the Newton equation. From the ensemble of final points we estimate $p(x)$ and $p(v)$ densities. Finally, the scenario (iii) assumes randomized motion with a fixed initial condition. In (iii) we have a fixed initial condition but the velocity is randomized: at random time instants $v(t_{i})\to-v(t_{i})$. For (iii), we see transient fingerprints of the initial condition which are sensitive to the waiting time distribution $p(\tau)$. At the same time, asymptotics $p(x)$ and $p(v)$ densities are robust with respect to the initial conditions on the constant energy curve and the waiting time density $p(\tau)$, but they depend on the total energy $E$.

Scenarios (i) – (iii) produce the same $p(x)$ and $p(v)$ distributions under the condition that the initial energy is the same, the observation time is long enough to reach the asymptotic regime. Exemplary Fig. 1 demonstrates the perfect agreement between results of numerical simulations of the scenario (iii) for $V(x)=x^{6}/6$ with stationary densities $p(x)$ and $p(v)$ given by Eqs. (5) and (6) which are derived using scenario (i). Therefore, in addition to decaying deterministic peaks, the main consequence of such randomization is the equivalence of ensemble and time averaging of motions starting with the fixed initial condition. This equivalence allows one to construct the exact distributions, see Eqs. (5) and (6).

Deviations from the obtained probability distributions can be easily produced by the addition of noise to Eq. (1). Presence of white Gaussian, Ornstein-Uhlenbeck or Markovian dichotomous noise results in the pumping of the energy to the system [29, 30, 31]. Consequently, the system’s trajectories in the phase space are dispersed and do not stay on closed orbits. Nevertheless, the trajectory $x(t)$ crosses $x=0$ infinitely many times [32]. The addition of noise without the introduction of the counterbalancing dissipative term introduces significant disturbances to the model. The alternative, more gentle scenario is to consider random initial conditions. In such a case, contrary to the fixed initial condition, the motion is not restricted to the single orbit, but a set of orbits determined by random initial conditions. After fixing the initial energy, the motion with hard velocity reversals is continued along the constant energy curve as in the scenario (iii). For random initial conditions, stationary densities can be calculated by the composition of probability densities

[TABLE]

and

[TABLE]

where $\lambda$ indicates a randomized parameter in $p(x)$ and $p(v)$ densities, while $p(\lambda)$ is the probability density of the randomized parameter. For simple distributions $p(\lambda)$ resulting densities given by Eqs. (7) and (8) can be calculated exactly.

Already for the uniform $p(\lambda)$ resulting formulas are complicated. Consequently, we present only a few representative figures and refer the reader to B for the full expressions. Figures 2 and 3 present sample $p(x)$ and $p(v)$ distributions for random initial conditions. In Fig. 2 the initial velocity is uniformly distributed over $[0.2,1.1]$ (left column) or $[0.9,1.1]$ (right column). Various rows correspond to different potentials $V(x)=x^{2}$ (top row) and $V(x)=x^{6}$ (bottom row). In Fig. 3, the potential is fixed to $V(x)=x^{4}$ and the randomized parameter is the initial total energy $E$ which is uniformly distributed over $[0.2^{2}/2,1.1^{2}/2]$ (top panel) $[0.9^{2}/2,1.1^{2}/2]$ (bottom panel). The selected support of energy distributions corresponds to minimal and maximal velocities considered in Fig. 2. Stationary densities corresponding to random initial conditions differ from the arcsine density corresponding to the fixed initial condition, see Fig. 1. For more complicated distributions $p(\lambda)$ stationary densities have to be calculated numerically.

3.2 2D deterministic and randomized motions

3.2.1 Deterministic motion

The Newton equation considered in Sec. 3.1 can be readily generalized to a higher number of dimensions

[TABLE]

Here we are interested in the case where $V(r)$ is a 2D, static single-well potential. As exemplary potentials, we use $V(r)=r^{n}/n$ with $n\in\{2,4\}$. Trajectories of Eq. (9) are bounded, i.e., they are located within a circle of finite radius $r_{\mathrm{max}}$. More precisely, they are located within the annulus defined by $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$ radii. As in 1D, the system is conservative, i.e., its energy is conserved, and the motion is continued along the constant energy curve. Moreover, in the absence of velocity reversals also the angular momentum $L$ is constant. According to the Bertrand theorem [33], Eq. (9) has closed orbits for $n=2$ (2D harmonic oscillator) and $n=-1$ (Kepler/Coulomb problem, not considered here) [34]. In the case of $n=2$, the marginal probability distributions can be found analytically. The marginal densities in the Cartesian coordinates, are given by

[TABLE]

where $A_{x}$ is the amplitude of motion along the $x$-axis. Using the transformation of variables and Eq. (30) one gets

[TABLE]

Formulas for $p_{y}(y)$ and $p_{v}(v_{y})$ are given by Eqs. (10) and (11) with $x$ replaced by $y$. The derivation of Eqs. (10) and (11) is included in C. Analogously, marginal densities can also be found in the spherical coordinates; appropriate formulas are included in D. Contrary to $n=2$, for $n=4$, solutions of Eq. (9) do not produce closed, non-circular, orbits. Moreover, there is no exact solution for $\vec{r}(t)$. Therefore, the approach used in the case of $n=2$ in appendices needs to be modified. In all cases, full and marginal densities can be easily obtained numerically, and some exact results can be constructed even for $n\neq 2$.

For $n=2$ the motion is performed along an ellipse. Therefore, the probability distribution $p(x,y)$ is defined on this ellipse only. The same positions are attained at equidistant times, which for the selected parameters differ by $T=2\pi$. The $p(x,v)$ density is shown in the top panel of Fig. 4. In Fig. 4 we use $r(0)=2,\dot{r}(0)=2.1,\varphi(0)=0$ and $\dot{\varphi}(0)=5.1$. The points of maximal probability are the same as predicted by the 2nd Kepler’s law. The probability density $p(x,y)$ is maximal in places where the velocity $v(x,y)$ is minimal, i.e., in the points which are farthest from the origin (center of the ellipse). The lower panels of Fig. 4 depict marginal densities $p_{x}(x)$ and $p_{x}(y)$ (middle row), $p(r)$ and $p(\varphi)$ (bottom row). Solid lines present theoretical formulas for $n=2$, see Eqs. (10), (D) and (44), which match results of computer simulations. The fingerprints of the 2nd Kepler’s law are also visible in $p(r)$ and $p(\varphi)$ distributions. The maxima of $p(\varphi)$ are located at points where $\dot{\varphi}$ is minimal. The angular velocity $\dot{\varphi}$ assumes the smallest values at the points which are the farthest from the origin (apocenter), i.e., at points corresponding to the maximal value of $r(t)$. By analogy: minimum of $p(\varphi)$ is located at the point where angular velocity is maximal, i.e., at the point which is the closest (pericenter) to the origin. Using the property $p(\varphi)\propto 1/\dot{\varphi}$, see Eq. (D), the distribution $p(\varphi)$ is constructed in D. The numerically estimated $p(\varphi)$ perfectly follows the theoretical predictions, see the right bottom panel of Fig. 4. Using Kepler’s law, it is also possible to predict the location of extrema of $p(r)$. Maxima of $p(r)$ are in places where $\dot{r}$ is minimal because $p(r)\propto 1/\dot{r}$, see Eq. (40). The radial velocity $\dot{r}$ is minimal ($\dot{r}=0$) at points which are farthest (apocenter) and closest (pericenter) to the origin. The $p(r)$ distribution is not symmetric, because at the apocenter the radial velocity changes faster than at the pericenter. Consequently, for both distributions $p(\varphi)$ and $p(r)$ most and least probable values of $r$ and $\varphi$ are defined by geometric properties of the orbit, i.e., by its apocenter and pericenter. At the same time, $p_{x}(x)$ and $p_{y}(y)$ distributions are of the arcsine type, see Eq. (10) and C. These distributions are symmetric along the $x=0$ ($y=0$) line. Here, again results of computer simulations perfectly follow theoretical densities, see the middle panel of Fig. 4.

The results for the quartic oscillator $V(r)=r^{4}/4$ with the same initial conditions as for $n=2$ in Fig. 4 are presented in Fig. 5. Top row presents a short, 2D trajectory and 2D $p(x,y)$ density while lower panels depict marginal densities $p_{x}(x)$ and $p_{y}(y)$ (middle row), $p(r)$ and $p(\varphi)$ (bottom row). Solid lines in the bottom row present theoretical formulas, see Eqs. (40) and (44), which confirm the computer simulations. Marginal densities $p_{x}(x)$ and $p_{y}(y)$, see the middle row, are not only the same but they are symmetric along $x=0$ ($y=0$) line.

For $n=4$, although the discussion of location of maximal values of $p(r)$ and $p(\varphi)$ distributions is the same as for $n=2$, the problem becomes more complex. On the one hand, Kepler’s law still holds. On the other hand, the orbit is no longer closed, see the left top panel of Fig. 5. Distances from the origin to apocenters ($r_{\mathrm{min}}$) and pericenters ($r_{\mathrm{max}}$) are fixed and well-defined, but they regularly change their location on appropriate circles. Furthermore, the time interval between consecutive visits to the apocenter ($r_{\mathrm{min}}$) and the pericenter ($r_{\mathrm{min}}$) is not fixed, but varies. Consequently, it is not possible to construct $p(r)$ and $p(\varphi)$ densities in an analogous way like for $n=2$. Nevertheless, $p(r)$ can be calculated using more general approach, see D where further differences between $p(r)$ and $p(\varphi)$ distributions are discussed. The most general formula for $p(r)$, see Eq. (44), reads

[TABLE]

where $U_{\mathrm{eff}}(r)$ is the effective potential

[TABLE]

and $L$ is the angular momentum. The system’s total energy $E$ and the angular momentum $L$ are determined by the initial condition, see D. Solid lines in left bottom panels of Figs. 4 and 5 correspond to Eq. (12). It is also possible to derive analytically the $p(\dot{\varphi})$ distribution, see Eq. (D). In this case, results of computer simulations nicely follow Eq. (D) (results not shown).

The top right panel of Fig. 5 presents the 2D histogram for the motion in the quartic potential. The presented histogram is constructed from the long yet finite trajectory. Within the histogram fingerprints of the deterministic trajectory are visible. Nevertheless, already a finite realization allows for perfect reconstruction of the marginal $p(r)$ and $p(\varphi)$ densities. With the increasing trajectory length, the internal part of the histogram becomes filled, but will not become uniform. If $p(x,y)$ was uniform on the disk $r_{\mathrm{min}}\leqslant r\leqslant r_{\mathrm{max}}$ the $p(r)$ density would be linear. This is definitely not the case, see the left bottom panel of Fig. 5 and Eq. (12). Therefore, even in the $t\to\infty$ limit $p(x,y)$ will not be uniform.

3.2.2 Randomized motion

Similarly to the 1D case, the motion defined by Eq. (9) is performed along the constant energy curve which for $n\neq 2$ (in general) is not closed [33]. Hard velocity reversals not only correspond to jumps in the phase space but also to change in the sign of the angular momentum $L$. Subsequent Figs. 6 and 7 depict main results of the analysis of the deterministic motion in 2D single-well potential of $r^{2}/2$ and $r^{4}/4$ type accompanied with hard velocity reversals. Velocity reversals are performed at random times which are the sum of random increments $\tau_{i}$ distributed according to a heavy-tailed probability density (the same protocol as in 1D). Analogously like in Dybiec et al. [5], we assume that $\tau_{i}=|\zeta_{i}|$, where $\zeta_{i}$ are independent identically distributed random variables following symmetric $\alpha$-stable density with the scale parameter equal to unity and the characteristic function given by

[TABLE]

The stability index $\alpha$ ($0<\alpha\leqslant 2$) describes power-law asymptotics of symmetric $\alpha$-stable densities. For large values of argument these densities decay as $|\zeta|^{-(\alpha+1)}$. Consequently, for $0<\alpha<1$ they are characterized by the diverging mean, while for $1<\alpha<2$ the variance diverges. Finally, in the limit of $\alpha=2$, they become normal distributions with all moments finite.

The 2D extension of the 1D model studied in Dybiec et al. [5] displays similar properties to the 1D model. For finite time and $\alpha$ small enough, $p_{x}(x)$, $p_{y}(y)$, $p(r)$ and $p(v)$ densities do not attain their stationary limits, which analogously like in 1D are robust with respect to the type of $p(\tau)$. This is clearly visible for $\alpha=0.1$ and $t=1000$ where violet squares in Figs. 6 and 7 are placed significantly below other curves because a substantial part of the probability mass is located in the peak corresponding to the fully deterministic dynamics performed without velocity reversals. Also for $\alpha=1$ (green bullets) some peaks are visible. Finally, for $\alpha>1$ the mean value $\langle\tau\rangle$ is finite and densities attain their stationary shapes, which are the same as for the fully deterministic motions studied in Figs. 4 and 5.

4 Summary and conclusions

The deterministic, frictionless, Newtonian motion in single-well 1D or 2D potentials is a universally studied example of a conservative system. Both in 1D and 2D systems the total energy is conserved, whereas in 2D the angular momentum is also constant. Very long observation of such a motion allows one to estimate the probability of recording a given velocity or position.

The deterministic motion can be disturbed in numerous ways; usually instantaneously breaking both momentum and energy conservation. In contrast, the studied here types of randomization preserve the energy conservation but break the angular momentum conservation. Example of such randomization is the random velocity reversal process [5]. Despite fixed initial conditions due to velocity reversals the averaging over a single long trajectory can be used interchangeably with the ensemble averaging. Accordingly, in 1D, it was possible to significantly extend results of Ref. [5]. Using more general, fully analytical approach, we have obtained general analytical formulas for $p(x)$ and $p(v)$ densities, which are valid for all single-well potential of $|x|^{n}$ type. In 1D stationary densities $p(x)$ and $p(v)$ are symmetric and $u$-shaped. The deviations from general $u$-shaped densities can be introduced by the additional randomization of the initial conditions, such randomizations are also traceable analytically (for some simple cases).

Finally, the 1D model has been extended to 2D setups. In contrast to 1D single-well potentials, 2D single-well potentials result in non-closed orbits (except for $n=2$) and quasi-periodic motion. Nevertheless, even in such a case, it is possible to derive exact formulas, at least for marginal $p_{x}(x),p_{y}(y),p_{v}(\dot{x}),p_{v}(\dot{y})$ and $p(r)$ densities which asymptotically are insensitive to velocity reversals.

This research was supported in part by PLGrid Infrastructure. Computer simulations have been performed at the Academic Computer Center Cyfronet, AGH University of Science and Technology (Kraków, Poland) under CPU grant “DynStoch”. Discussions with Ralf Metzler and Aleksei Chechkin are greatly acknowledged.

Appendix A Derivation of formula for $p(x)$ and $p(v)$

In Dybiec et. al. [5] formulas for $p(x)$ and $p(v)$ have been derived using semi-analytical and numerical methods for the potential of $|x|^{n}$ type with $n\in\{1,2,4,\infty\}$. Here, we present general, fully analytical, derivation that is valid for any exponent $n$.

For the 1D periodic motion the probability of finding a particle at the point $x$ is proportional to the time spent in the vicinity of this point

[TABLE]

where $T$ is the period of the motion [28] in $V(x)=\kappa|x|^{n}$ ($\kappa>0,n>0$)

[TABLE]

The velocity $v=\dot{x}$ can be calculated from the energy conservation principle

[TABLE]

resulting in

[TABLE]

The total energy $E$ is determined by the initial condition $E=\frac{1}{2}mv_{0}^{2}+\kappa|x_{0}|^{n}$. Finally, the formula for the $p(x)$ density reads

[TABLE]

The particle can be most easily found at points where it spends a long time, see Eq. (15), i.e., at points where the velocity is minimal. These points correspond to places at which the velocity is softly reversed, i.e., to the points which are in the maximal distance from the origin, see Eq. (19).

The $p(v)$ density can be obtained by the transformation of variables

[TABLE]

From the energy conservation

[TABLE]

and

[TABLE]

Finally, after substituting Eq. (22) in Eq. (20) and using the fact that $E-\kappa|x|^{n}=\frac{1}{2}mv^{2}$ one obtains

[TABLE]

Maxima of the $p(v)$ density are located at points where the acceleration is minimal, i.e., the velocity is maximal. At the origin ($x=0$), the velocity is maximal while the acceleration is minimal. Consequently, at the origin, the $p(x)$ is minimal, while $p(v)$ is maximal.

Appendix B Uniform energy distribution

In simplest cases formulas (7) and (8) can be integrated analytically. For instance, if the randomized parameter is the total energy $E$ and it is uniformly distributed over $[E_{1},E_{2}]$ one can calculate $\tilde{p}(E,x)=\int p(x|E)p(E)dE$ and $\tilde{p}(E,v)=\int p(v|E)p(E)dE$ resulting in

[TABLE]

for $E-\kappa\left|x\right|^{n}\geqslant 0$, otherwise $\tilde{p}(E,x)=0$. Analogously

[TABLE]

for $E-\frac{m}{2}v^{2}\geqslant 0$, otherwise $\tilde{p}(E,v)=0$.

Finally, the stationary distributions $p(x)$ and $p(v)$ read

[TABLE]

and

[TABLE]

Appendix C Motion in $V(r)=\kappa r^{2}$

For $V(r)=\kappa r^{2}$ ($\kappa>0$) orbits of the motion in 2D are closed. Moreover, in the Cartesian coordinates, Eq. (9) separates into two independent equations

[TABLE]

which can be easily solved. For instance, the formula for $x(t)$ reads

[TABLE]

where $\omega^{2}=2\kappa/m$, $\Delta_{x}=-\arctan{\frac{\dot{x}_{0}}{\omega x_{0}}}$ and $A_{x}=\sqrt{x_{0}^{2}+\dot{x}^{2}_{0}/\omega^{2}}$. The analogous formula holds for $y(t)$. The marginal density $p_{x}(x)$ can be found in an analogous way like $p(x)$ in A

[TABLE]

The velocity $\dot{x}$ as a function of the position $x$ can be calculated from the condition

[TABLE]

which allows to calculate

[TABLE]

Using the transformation of variables and Eq. (30) one gets

[TABLE]

Marginal densities $p_{y}(y)$ and $p_{v}(\dot{y}=v)$ are given by Eqs. (31) and (32) with $x$ replaced by $y$, i.e., $A_{x}\to A_{y}$, $x\to y$ and $\dot{x}\to\dot{y}$. These densities are arcsine-like and even.

Appendix D 2D Motion in $V(r)=\kappa r^{n}$

The equation of motion (9) in the central potential $V(r)=\kappa r^{n}$ ($\kappa>0$ and $n>0$) can be rewritten in the polar coordinates $(r,\varphi)$ as

[TABLE]

where $U_{\mathrm{eff}}(r)$ is the effective potential

[TABLE]

and $L$ is the angular momentum

[TABLE]

The angular momentum is conserved and its value is determined by the initial condition

[TABLE]

For $n=2$ and the special case of $\omega^{2}=2\kappa/m=1$, $m=1$ the solution of Eq. (33) can be found analytically

[TABLE]

From Eq. (37) and the conservation of angular momentum $L$ the angular velocity $\dot{\varphi}(t)$ together with $\varphi(t)$ and its inversion $\varphi^{-1}(t)$ can be evaluated. Here, for brevity, we present only $\varphi(t)$

[TABLE]

Note, that initially $\varphi(0)=\arctan{\left[\frac{mr_{0}\dot{r}_{0}}{L}\right]}$. For such an initial condition, the probability distribution $p(\varphi)$ can be obtained analytically

[TABLE]

Analogously to $p(\varphi)$, the marginal density $p(r)$ can be found from

[TABLE]

In Eq. (40) there is the $\propto$ sign. The origin of this proportionality is visible for $n=2$. It originates due to the fact that the same radii $r$ are taken more than once during one period of motion. More precisely, minimal and maximal values of $r$ are recorded twice, see below, while intermediate values are taken four times. The formula for $p(r)$ can be derived with the help of Eq. (40) and Eq. (37). Unfortunately, it is not the most straightforward approach and the constructed formula is quite lengthy. The more convenient approach is to utilize the energy conservation and the angular momentum conservation, see below.

For $n\neq 2$ orbits do not close; nevertheless, the motion is bounded. The motion is restricted to $r(t)\in[r_{\mathrm{min}},r_{\mathrm{max}}]$ where $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$ are determined by the condition

[TABLE]

where $E$ is the total (conserved) energy and $U_{\mathrm{eff}}(r)$ is the effective potential given by Eq. (34). Due to boundness of the motion marginal densities $p_{x}(x)$ and $p_{y}(y)$ as well as $p_{v}(v_{x})$ and $p_{v}(v_{y})$ have the same support and functional dependence. Please note however, that for $n\neq 2$, the marginal density $p(\varphi)$ cannot be constructed using the approach applied in Eq. (D). First of all, orbits are not closed; thus the period of the motion is not defined. After increasing $\varphi(t)$ by $2\pi$ the radius $r(t)$ takes a different value, more precisely $r(\varphi)\neq r(\varphi+2\pi)$. Therefore, due to the conservation of the angular momentum

[TABLE]

the mapping $\varphi\to\dot{\varphi}$ is not unique. Nevertheless, the $p(\varphi)$ density can be estimated numerically. From computer simulations it seems that $p(\varphi)$ is uniform. At the same time the marginal density $p(r)$ still can be constructed with the help of Eq. (40). Despite lack of periodicity, due to the energy and momentum conservations

[TABLE]

the mapping $r\to\dot{r}$ is unique and Eq. (40) is valid. The formula for $p(r)$ reads

[TABLE]

where $U_{\mathrm{eff}}(r)$ is given by Eq. (34) and $E$ and $L$ are determined by the initial condition. Eq. (44) gives the general formula for the $p(r)$ density for any 2D single-well potential. Practically, the only difficulty is the calculation of the $r_{\mathrm{min}}$, $r_{\mathrm{max}}$ and the normalization constant $N$

[TABLE]

Furthermore, using the transformation of variables and the relation between the angular momentum and the radial velocity from Eq. (44) it is possible to calculate

[TABLE]

The support of $p(\dot{\varphi})$ density is determined by the $L=mr^{2}\dot{\varphi}$ condition. In accordance with Kepler’s law, the minimal (maximal) radial velocity corresponds to the maximal — $r_{\mathrm{max}}$ (minimal — $r_{\mathrm{min}}$) radius.