*Article* **The Relativistic Harmonic Oscillator in a Uniform Gravitational Field**

**Michael M. Tung**

Instituto Universitario de Matemática Multidisciplinar, Universitat Politècnica de Valencia, Camino de Vera, s/n, 46022 Valencia, Spain; mtung@mat.upv.es

**Abstract:** We present the relativistic generalization of the classical harmonic oscillator suspended within a uniform gravitational field measured by an observer in a laboratory in which the suspension point of the spring is fixed. The starting point of this analysis is a variational approach based on the Euler–Lagrange formalism. Due to the conceptual differences of mass in the framework of special relativity compared with the classical model, the correct treatment of the relativistic gravitational potential requires special attention. It is proved that the corresponding relativistic equation of motion has unique periodic solutions. Some approximate analytical results including the next-to-leadingorder term in the non-relativistic limit are also examined. The discussion is rounded up with a numerical simulation of the full relativistic results in the case of a strong gravity field. Finally, the dynamics of the model is further explored by investigating phase space and its quantitative relativistic features.

**Keywords:** relativistic harmonic oscillator; kinematics of a particle; special relativity; nonlinear problems in mechanics; equations of motion in gravitational theory

**MSC:** 70B05; 83A05; 70K42; 70K99; 83C10

### **1. Introduction**

The harmonic oscillator or mechanical spring, implementing Hooke's Law, is one of the standard textbook examples for introducing the student to Newtonian mechanics [1]. Treating classical motion under a Hookean potential is simplest, in spite of additional difficulties when, e.g., velocity-dependent forces (as friction) are added. Remarkably, the relativistic counterpart of an oscillator or a pendulum—which approximates to a harmonic oscillator for small amplitudes—stationed within some supplementary force field has so far been dealt with only scarcely.

A detailed discussion of relativistic effects on a simple pendulum without any additional forces has been carried out by Erkal in 2000 [2]. In 2008, Torres shows that the relativistic pendulum with friction possesses periodic solutions which are absent in the classical case [3]. In a more recent publication of this area of research, in 2017, de la Fuente and Torres focuses on relativistic extensions for the motion of the harmonic oscillator from the view of the oscillating body, but without including any gravitational effects [4]. In all these publications an appropriate laboratory frame is chosen where characteristic, preferred points are fixed, i.e., the suspension points of the pendulum and the spring representing the oscillator, respectively.

The present work fills an outstanding gap in the existing research literature by examining the relativistic effects of a harmonic oscillator in a uniform gravitational field, adopting and extending the approach by Goldstein and Bender [5]. In the classical model, the maximum velocity of the mass in motion can be arbitrarily large depending on the displacement with respect to the equilibrium point of the spring. Relativistic mechanics will adjust this behavior by only allowing a maximum speed less than the speed of light, as accelerating a mass to higher velocities will in like manner increase its inertial mass

**Citation:** Tung, M.M. The Relativistic Harmonic Oscillator in a Uniform Gravitational Field. *Mathematics* **2021**, *9*, 294. https://doi.org/10.3390/ math9040294

Academic Editors: Rami Ahmad El-Nabulsi; Andrea Scapellato and Lucas Jódar Received: 9 January 2021 Accepted: 1 February 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(equivalent to an increase in kinetic energy). Although the oscillating mass point is not considered to generate its own gravitational field and the uniform gravitational field shall be caused by an external object, any variable inertial mass will also be subject to gravity. This will inescapably lead to further nontrivial complications when studying relativistic effects of gravity for the harmonic oscillator—even from the perspective of special relativity without considering general relativity. Nonetheless, the present model may provide a legitimate and satisfactory approximation for an oscillator close to a black hole, as on a small scale the Schwarzschild spacetime at a sufficiently large distance from the event horizon (where gravitational tidal effects can be ignored) approximates well a strong uniform gravitational field.

The core of this paper is organized as follows. In Section 2, after a short introduction to the variational principle [6,7] for the classical case of the harmonic oscillator, the Euler– Lagrange formalism in the framework of special relativity [8,9] is used to set up the equations of motion for a harmonic oscillator which will be subject to an external uniform gravitational field and measured by an observer in a laboratory where the suspension point of the spring is fixed. Generalizing from the classical to the relativistic regime is not as trivial as it appears at first sight due to the particular, distinct nature of relativistic mass—mass which attains a dynamical quality—and the fact that variable kinetic energy itself is equivalent to additional mass which is susceptible to the external gravitational field. Special care has to be taken to take these entirely relativistic effects into account.

Therefore, Section 3 concentrates on the full derivation of the correct relativistic potential for the uniform gravitational field surrounding the harmonic oscillator. As approximation in the case of weak gravitational fields, we consider the Taylor expansion of the potential in the non-relativistic limit and some of its particular properties. Furthermore, we examine the full relativistic results for the potential with strong gravity and, in particular, identify its physically allowed regions. Although we are able to derive the relativistic gravitational potential for the case at hand, and in closed analytical form, the final results for the equation of motion become intangible for analytical evaluation.

In Section 4, we perform the numerical integration of the equation of motion to simulate the dynamics of this model and explore some significant characteristics of the system. In general, we prove that the equation of motion for the relativistic harmonic oscillator also has unique periodic solutions. In the strong gravity case, we compare relativistic with classical estimates for the oscillating amplitudes. For further analysis, we present the corresponding phase-space trajectories and discuss its most prominent characteristics.

#### **2. Variational Principle and Equation of Motion**

Robert Hooke (1635–1703) first pointed out that the mathematical description for small oscillations of a body with mass *m*<sup>0</sup> > 0 attached to an elastic spring with position *<sup>x</sup>* = *<sup>x</sup>*(*t*) takes the form: *<sup>m</sup>*0*x*¨ = <sup>−</sup>*k x*. The positive constant *<sup>k</sup>* <sup>&</sup>gt; 0 depends on the elastic properties of the spring in question. As a natural length scale serves the length of the spring at its maximum elongation, denoted as > 0. This mechanical system is termed "harmonic oscillator". Such systems are of utmost relevance in physics and engineering, as any mass particle subject to a force in stable equilibrium will effectively operate as a harmonic oscillator for small fluctuations—small fluctuations being displacements with only a fraction of length . Additional importance emerges in the dynamics of a continuous classical field as it may be formulated as the dynamics of an infinite number of harmonic oscillators. Furthermore, the quantum harmonic oscillator describes some of the most important model systems in quantum mechanics.

Already for the elementary classical case of a spring extended in a uniform Newtonian gravitational field, the most efficient and powerful approach is the framework of

Lagrangian mechanics [6,7], which is based on a variational principle. Accordingly, the deterministic equations of motion will result from the following principle of least action

$$\delta \int \text{dt } \mathbf{L} = \delta \int \text{dt} \left[ \frac{1}{2} m\_0 \dot{\mathbf{x}}^2 - \frac{1}{2} k \mathbf{x}^2 - m\_0 \mathbf{g} \, \mathbf{x} \right] = \mathbf{0} \tag{1}$$

by varying over all possible paths *<sup>x</sup>* = *<sup>x</sup>*(*t*) and keeping the end points fixed. From Equation (1), it becomes clear that we exclude the oscillating mass point (which has negligible mass) as possible source of a gravitational field. For practical purposes and as is customary, the position *x* of the mass point and the external uniform gravitational field with strength *g* > 0 are measured in the laboratory frame, where the suspension of the spring at one end is fixed.

The integrand in Equation (1) is the Lagrangian function, *L*, which includes kinetic energy *<sup>T</sup>* = <sup>1</sup> <sup>2</sup>*m*0*x*˙ 2, the spring potential *Vs*(*x*) = <sup>1</sup> <sup>2</sup> *kx*2, and the gravitational potential *Vg*(*x*) = *<sup>m</sup>*0*gx*. Note that we consider spring elongations with respect to position *<sup>x</sup>* = <sup>0</sup> which is also the suspension point of the spring. In this laboratory frame, with the suspension point at rest, the oscillating device is stationed within a uniform gravitational field (determined by the gravitational constant *g* > 0) in such a way that both are aligned. This dynamical system is one-dimensional, and the corresponding Euler–Lagrange equation amounts to solving a simple linear second-order differential equation deriving from

$$
\left(\frac{d}{dt}\frac{\partial}{\partial \dot{x}} - \frac{\partial}{\partial x}\right)L = 0.\tag{2}
$$

By substituting the Lagrangian *L* from Equation (1) into Equation (2), it is straightforward to reproduce the well-known general solution in closed analytical form:

$$\mathbf{x}(t) = \mathbf{C}\_1 \cos\left(\sqrt{\frac{k}{m\_0}} \, t\right) + \mathbf{C}\_2 \sin\left(\sqrt{\frac{k}{m\_0}} \, t\right) - \frac{m\_0 \mathbf{g}}{k},\tag{3}$$

where *C*<sup>1</sup> and *C*<sup>2</sup> are the two integration constants depending on the initial values for the differential equation.

However, the classical result, Equation (3), does not contemplate strong gravitational fields and when velocities *x*˙ get closer to the speed of light *c* > 0. Furthermore, in special relativity the mass is a dynamical quantity, dependent on the relative velocity *x*˙ of the observer, such that *<sup>m</sup>* = *<sup>γ</sup>m*0, where the usual relativistic factor is *<sup>γ</sup>*(*x*˙) = 1/ <sup>√</sup><sup>1</sup> <sup>−</sup> *<sup>x</sup>*˙2/*c*2. Consequently, Equation (3) will utterly fail in giving a faithful description of the physical effects in the relativistic domain.

In order to generalize to a correct description in the relativistic domain, the best starting point is to modify the classical principle of least action, Equation (1). As expected, the rest mass *m*<sup>0</sup> in Equation (1) will have to be divided by the factor *γ* to correctly incorporate both rest mass and kinetic energy. Moreover, the spring potential *Vs* is unaltered, and the relativistic gravitational potential *Vg*, however, is hitherto undetermined. Therefore, we postulate the relativistic Lagrangian

$$L(\mathbf{x}, \dot{\mathbf{x}}) = -\frac{m\_0 c^2}{\gamma(\dot{\mathbf{x}})} - \frac{1}{2} k \mathbf{x}^2 - V\_{\mathcal{S}}(\mathbf{x}),\tag{4}$$

which readily yields the relativistic action principle

$$\delta \int dt \left[ -m\_0 c^2 \sqrt{1 - \frac{\dot{\mathbf{x}}^2}{c^2}} - \frac{1}{2} k \mathbf{x}^2 - V\_\delta(\mathbf{x}) \right] = 0. \tag{5}$$

Note that *Vg*(*x*), as stressed before, gives the relativistic gravitational potential as measured in the laboratory frame with fixed strength *g* > 0 and as a function of varying spring elongation *x*.

Complications for identifying *Vg* in Equation (5) arise, because mass in special relativity will change with its variable kinetic energy and thus simultaneously cause a change of the mass which is subject to the external gravitational field. This effect can be taken into account by using the relativistic ansatz

$$W\_{\mathcal{J}}(\mathbf{x}) = \int d\mathbf{x} \, m \, \mathbf{g} = m\_0 \mathbf{g} \int \frac{d\mathbf{x}}{\sqrt{1 - \frac{\mathbf{g}^2}{c^2}}}. \tag{6}$$

For a rigorous general definition of this potential see Goldstein & Bender [5].

As *Vg*, and thus *L*, do not explicitly depend on time, according to Noether's theorem total energy must be conserved. For this purpose, we carry out the Legendre transformation [6,7] of Equation (4) yielding the Lagrangian energy function

$$E(\mathbf{x}, \dot{\mathbf{x}}) = \frac{m\_0 c^2}{\sqrt{1 - \frac{\mathbf{x}^2}{c^2}}} + \frac{1}{2} k \mathbf{x}^2 + V\_{\mathcal{S}}(\mathbf{x}) =: E,\tag{7}$$

which is just the constant total energy *E* of the system. Without loss of generality, but for convenience, we assume for the remainder of the derivation that at position *<sup>x</sup>* = 0 the mass particle be at rest, i.e., *<sup>E</sup>* = *<sup>m</sup>*0*c*2, or equivalently, we will assume that the following initial conditions hold:

$$
\dot{\mathbf{x}}(0) = 0 \quad \text{and} \quad \dot{\mathbf{x}}(0) = 0. \tag{8}
$$

Furthermore, this immediately implies via Equation (7) that

$$V\_{\S}(0) = 0.\tag{9}$$

Observe that condition Equation (9) is chosen to agree with the definition of the classical potential and is nothing more than just fixing the arbitrary, and unphysical, integration constant in Equation (6).

Applying the Euler–Lagrange formalism to Equation (5) produces the following relativistic equation of motion,

$$\frac{d}{dt}\left(\frac{m\_0 \dot{\mathbf{x}}}{\sqrt{1 - \frac{\dot{\mathbf{x}}^2}{c^2}}}\right) + k\mathbf{x} + \frac{d}{d\mathbf{x}}V\_{\mathcal{S}}(\mathbf{x}) = \mathbf{0},\tag{10}$$

where *Vg*(*x*) is still unknown and needs to be determined. The next section focuses on uncovering the explicit relativistic form of gravitational potential *Vg*(*x*).

#### **3. Relativistic Potential for Uniform Gravitational Field**

As already stressed, the problem of dealing with the relativistic model of the harmonic oscillator in a uniform gravitational field is appreciably more complicated than the classical case. The core problem originates from the fundamentally different concept of mass in the classical or the relativistic description of the physical phenomena. The ansatz for the gravitational potential given in Equation (6) naturally considers relativistic corrections of a mass in relative motion with respect to the laboratory.

Moreover, from energy conservation, viz. Equation (7), with *<sup>E</sup>* = *<sup>m</sup>*0*c*2, and from solving for the relativistic factor *γ*, it directly follows that

$$\gamma = \frac{1}{\sqrt{1 - \frac{\chi^2}{c^2}}} = \frac{m\_0 c^2 - \frac{1}{2} k x^2 - V\_\%(x)}{m\_0 c^2} \ge 1,\tag{11}$$

which constrains the physically admissible range of potential *Vg*.

Now, substituting the relativistic factor of Equation (11) into Equation (6) gives the integral equation

$$V\_{\mathcal{S}}(\mathbf{x}) = \frac{\mathcal{S}}{c^2} \int d\mathbf{x} \left[ m\_0 c^2 - \frac{1}{2} k \mathbf{x}^2 - V\_{\mathcal{S}}(\mathbf{x}) \right],\tag{12}$$

or correspondingly, the differential equation

$$\frac{d}{d\mathbf{x}}V\_{\mathcal{S}}(\mathbf{x}) = \frac{\mathcal{S}}{c^2} \left[ m\_0 c^2 - \frac{1}{2} k \mathbf{x}^2 - V\_{\mathcal{S}}(\mathbf{x}) \right] > 0, \quad V\_{\mathcal{S}}(\mathbf{0}) = \mathbf{0},\tag{13}$$

where we have considered integration condition Equation (9). Furthermore, we indicated the correct sign of gradient *dVg*/*dx*, since the gravitational force, *Fg* <sup>=</sup> <sup>−</sup>*dVg*/*dx*, is oriented downwards the *x*-axis. It is important to notice that Equation (13) contains the constant *k* > 0, considering that obviously the elastic property of the spring affects the speed of the oscillating mass and enters *Vg* via the relativistic factor *γ*. This is in full agreement with previous published results using a similar approach, e.g., for the relativistic pendulum with a gravitational potential energy also depending on length of the pendulum, viz. (see in [2], Equation (6)).

Equation (13) is an ordinary differential equation of type *f* (*x*) = *<sup>a</sup>* + *bx*<sup>2</sup> + *c f*(*x*) which can easily be solved. The solution is

$$V\_{\mathcal{S}}(\mathbf{x}) = \left[m\_0 c^2 - k\left(\frac{c^2}{\mathcal{S}}\right)^2\right] \left(1 - e^{-\frac{\mathcal{S}}{c^2}\mathbf{x}}\right) + \frac{1}{2}k\left(\frac{2c^2}{\mathcal{S}} - \mathbf{x}\right)\mathbf{x} \tag{14}$$

with derivative

$$\frac{d}{dx}V\_{\mathcal{S}}(\mathbf{x}) = \left[m\_0 \mathbf{g} - k \frac{c^2}{\mathcal{S}}\right] e^{-\frac{\mathcal{S}}{c^2} \mathbf{x}} + k \left(\frac{c^2}{\mathcal{S}} - \mathbf{x}\right). \tag{15}$$

Some checks of Equation (15) are in order: Note that in the absence of any spring (*<sup>k</sup>* = 0), the result *Fg* <sup>=</sup> <sup>−</sup>*<sup>V</sup> <sup>g</sup>* <sup>=</sup> <sup>−</sup>*m*0*g e*<sup>−</sup> *<sup>g</sup> <sup>c</sup>*<sup>2</sup> *<sup>x</sup>* is recovered, which is in full agreement with the calculations by Goldstein and Bender [5]. As a consequence, the classical result, *Fg* <sup>=</sup> <sup>−</sup>*m*0*<sup>g</sup>* is obtained for *<sup>x</sup>* = 0, before the mass particle is set in motion, viz. Equation (8). Moreover, in the weak gravity domain (*g*/*c*<sup>2</sup> 1), to lowest order the gravitational force obviously has to be independent of spring constant *<sup>k</sup>*, and it is *Fg* ≈ −*m*0*<sup>g</sup>* <sup>+</sup> <sup>O</sup> *g*/*c*<sup>2</sup> , with being the natural length scale of the spring, viz. Section 2. However, as gravity becomes stronger, the harmonic force will entangle with gravity in the potential *Vg* due to the subtle relativistic effects already mentioned. To see this, we expand Equation (14) in a power series and obtain the expansion of *Vg* up to first order in the dimensionless scale parameter *g*/*c*2:

$$V\_{\mathcal{S}}(\mathbf{x}) = m\_0 \mathbf{g} \, \mathbf{x} - \frac{k \mathbf{x}^3}{6\ell} \left(\frac{\mathbf{g}\,\ell}{c^2}\right) + \mathcal{O}\left(\frac{\mathbf{g}\,\ell}{c^2}\right)^2. \tag{16}$$

Observe that the quantity *g*/*c*<sup>2</sup> is a Lorentz scalar, representing an invariant for all inertial frames. Obviously, the speed of light *c* is a Lorentz scalar, and *g* being an acceleration is measured the same in all inertial frames with relative motion to the rest/laboratory frame. Thus, the ratio *g*/*c*<sup>2</sup> is also invariant under Lorentz transformations.

In Equation (16), the term of order O *g*/*c*<sup>2</sup> already contains spring constant *k*. Therefore, in first approximation for weak gravity, the odd powers indicate a symmetric result, more precisely rotational symmetry with respect to rotations of 180◦ about the origin, see Figure 1. In the full relativistic regime, including all orders of the expansion, symmetry is broken. Observe also that in this expansion of the gravitational potential, the additional next-to-leading-order term represents the main correction of special relativity. It bears some similarity with the post-Newtonian approximation in general relativity. However, here in our approach—within the framework of special relativity—the underlying spacetime is of course Minkowskian, and thus flat.

**Figure 1.** In the weak gravity case, when *<sup>g</sup>*/*c*<sup>2</sup> 1, the approximation for the relativistic gravitational potential *Vg* in the non-relativistic limit is most suitable (blue curve), see Equation (16). This approximate result is symmetric with respect to rotations of 180◦ about the origin—a symmetry property which is lost in the full relativistic case. The approximate model is only physically valid well within the interval [−*x*0, *<sup>x</sup>*0], where *<sup>x</sup>*<sup>0</sup> <sup>=</sup> /2*m*<sup>0</sup> *<sup>k</sup> c*, so that sign-flips of the force field are avoided. The non-relativistic, classical case is indicated by the straight gray line.

Figure 2 displays the full relativistic result for *Vg* in a strong gravitational field. Here, to achieve *<sup>g</sup>*/*c*<sup>2</sup> = <sup>1</sup> <sup>m</sup>−<sup>1</sup> (measured in physical units m−1), all parameters are set to unity, except for the spring constant which we chose to be *<sup>k</sup>* = <sup>2</sup> kg/s2 (measured in physical units kg/s2). We also represent the gradient, *V <sup>g</sup>* <sup>=</sup> *dVg*/*dx*, for any position *<sup>x</sup>* <sup>∈</sup> [−1, 1]. Note that the gravitational force therefore will flip sign at *x*0, satisfying *V <sup>g</sup>*(*x*0) = 0 for Equation (13), and is given by

$$\mathbf{x}\_0 = \frac{c^2}{\mathcal{S}} \left[ \mathcal{W}\_0 \left( \frac{\frac{m\_0 \mathcal{S}^2}{kc^2} - 1}{c} \right) + 1 \right],\tag{17}$$

where *W*<sup>0</sup> is the principal branch of Lambert's *W* function (see in [10], §4.13). By including higher orders up to <sup>O</sup>(*g*/*c*2)<sup>3</sup> in Equation (16), a reasonably good approximation for Equation (17) is obtained: *<sup>x</sup>*<sup>0</sup> <sup>≈</sup> *<sup>m</sup>*<sup>0</sup> *<sup>g</sup> k* / 2*c*2*k <sup>m</sup>*<sup>0</sup> *<sup>g</sup>*<sup>2</sup> <sup>+</sup> <sup>1</sup> <sup>−</sup> <sup>1</sup> .

For the data in Figure 2, it is *x*<sup>0</sup> ≈ 0.77 m, and thus all estimates for *x* > *x*<sup>0</sup> are unphysical. However, with initial conditions Equation (8) the mass particle will move only on the negative axis, *x* ≤ 0, and thus will safely be in the physical region.

**Figure 2.** The relativistic gravitational potential, *Vg*(*x*), for an oscillating body posted in a strong uniform gravitational field. For the graphical representation we use *<sup>m</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> kg, *<sup>g</sup>* <sup>=</sup> <sup>1</sup> kg m/s2, and *<sup>k</sup>* = <sup>2</sup> kg/s2. Further, the speed of light is normalized to *<sup>c</sup>* = <sup>1</sup> m/s. The dashed line marks below the physically allowed region for *Vg*, viz. Equation (20). The gradient, *V <sup>g</sup>*(*x*), is also shown.

#### **4. Model Dynamics and Numerical Simulation**

With the exact analytical result for the relativistic gravitational potential *Vg*, given by Equation (14), we are now in the position to complete the description of the dynamics of the model at hand. Substituting its derivative *V <sup>g</sup>*, given by Equation (15), into the equation of motion, Equation (10), readily yields the Euler–Lagrange equation—a nonlinear secondorder differential equation—which governs the physical system

$$
\gamma^3 m\_0 \ddot{\mathbf{x}} + \left[ m\_0 \mathbf{g} - k \frac{c^2}{\mathcal{S}} \right] e^{-\frac{\mathcal{L}}{c^2} \mathbf{x}} + k \frac{c^2}{\mathcal{S}} = 0,\tag{18}
$$

where *γ* and *e* − *g <sup>c</sup>*<sup>2</sup> *<sup>x</sup>* can be eliminated via Equations (11) and (14), respectively. After some lengthy but straightforward simplification, we arrive at the following equivalent and for numerical implementation more convenient form

$$a(\mathbf{x})\ddot{\mathbf{x}} + \beta(\mathbf{x})\mathbf{g} = \mathbf{0},\tag{19a}$$

$$\begin{cases} a(\mathbf{x}) = \left(1 - \frac{\frac{1}{2}k\mathbf{x}^2 + V\_{\mathcal{S}}(\mathbf{x})}{m\_0 c^2}\right)^3, \\\\ \beta(\mathbf{x}) = \frac{\frac{1}{2}k\left(\frac{2\mathcal{L}^2}{\mathcal{S}} - \mathbf{x}\right)\mathbf{x} + m\_0 c^2 - V\_{\mathcal{S}}(\mathbf{x})}{m\_0 c^2},\end{cases}\tag{19b}$$

where *Vg* is given by Equation (14). Observe that *α*, *β* are dimensionless factors, and in particular it is *<sup>α</sup>*(0) = *<sup>β</sup>*(0) = 1, such that *<sup>x</sup>*¨ = <sup>−</sup>*<sup>g</sup>* at the initial position *<sup>x</sup>* = 0, as is expected. Furthermore, we obtain *<sup>x</sup>*¨ <sup>=</sup> <sup>−</sup>*<sup>g</sup>* at all positions, when *<sup>k</sup>* <sup>=</sup> 0 and *Vg* <sup>≡</sup> 0. This represents the non-relativistic case in the absence of harmonic forces, that is, classical free fall. Similarly, for *k* > 0 and *Vg* ≡ 0, it is easily checked that the result reduces to that of a classical spring in a gravitational field: *<sup>m</sup>*0*x*¨ = <sup>−</sup>*kx* <sup>−</sup> *<sup>m</sup>*0*<sup>g</sup>* with general solutions Equation (3). Anyhow, Equations (19a) and (19b) embraces the full relativistic case provided that constraint Equation (11) is satisfied. The constraint Equation (11) may simply be rewritten as

$$V\_{\mathcal{S}}(\mathbf{x}) \le -\frac{1}{2}k\mathbf{x}^2. \tag{20}$$

Therefore, the physically relevant gravitational potential *Vg* always has to lie below this concave down parabola. Figure 2 shows that this will approximately hold for the range −*x*<sup>0</sup> ≤ *x* ≤ 0, where no sign-flip for the gravitational force occurs. Recall that with the aforementioned parameters (all set to unity except for *<sup>k</sup>* = <sup>2</sup> kg/s2 ), we found *x*<sup>0</sup> ≈ 0.77 m, *viz.* Equation (17).

The particular structure of the relativistic equation of motion, Equations (19a) and (19b), also implies existence and uniqueness of periodic solutions. During the past two decades, considerable progress has been made in the study of second-order differential equations and the periodic properties of their solutions [11]. For a closer analysis, we recast Equation (19a) into the form

$$\ddot{x}(t) = -\frac{\beta(x)}{\kappa(x)}g =: f(x),\tag{21}$$

where in this case the continuous function *<sup>f</sup>* (assuming *<sup>α</sup>*(*x*) = 0) does not explicitly depend on the variables *<sup>t</sup>* or *<sup>x</sup>*˙. Moreover, recall that *<sup>f</sup>*(*x*(0)) = <sup>−</sup>*g*. To guarantee existence and uniqueness, we have to check that *f* is bounded by two continuous functions for all *<sup>t</sup>* <sup>∈</sup> [0, *<sup>T</sup>*], with period *<sup>T</sup>* <sup>&</sup>gt; 0 [11]. For this purpose, we also use the following general results for the limits of Equation (19b) and its derivatives:

$$\lim\_{x \to \infty} a = -\infty, \qquad \lim\_{x \to \infty} \beta = \frac{k}{m\_0} \left(\frac{c}{\mathcal{g}}\right)^2,$$

$$\lim\_{x \to \infty} a' = -\infty, \qquad \lim\_{x \to \infty} \beta' = 0,\tag{22}$$

$$\text{and} \qquad \lim\_{x \to \infty} \frac{a'}{a} = 0.$$

Now, we can conclude that

$$\lim\_{\alpha \to \infty} f' = \lim\_{\alpha \to \infty} \frac{\alpha' \beta - \alpha \beta'}{\alpha^2} g = 0. \tag{23}$$

As *f* (*x*) and *<sup>x</sup>*(*t*) are continuous, *<sup>f</sup>* necessarily has to be bounded by a lower and an upper continuous function for *t* ≥ 0, and we affirm that the equation of motion, Equations (19a) and (19b), has unique periodic solutions.

The series expansion in the non-relativistic limit of Equations (19a) and (19b) up to first order in the dimensionless parameter *g*/*c*<sup>2</sup> gives

$$2m\ddot{\mathbf{x}} + m\mathbf{q}\mathbf{g} + k\mathbf{x} + \left(\frac{2m\_0\mathbf{g}}{\ell}\mathbf{x} + \frac{4k}{\ell}\mathbf{x}^2 + \frac{3k^2}{2m\_0\mathbf{g}\ell}\mathbf{x}^3\right)\frac{\mathbf{g}\ell}{\mathbf{c}^2} = \mathbf{0},\tag{24}$$

and analytical integration is possible but will already yield a rather convoluted list of elliptic integrals. Therefore, for the full relativistic result—including all higher orders—numerical integration is the most appropriate, if not the only possible, open pathway.

For the actual numerical integration of Equations (19a) and (19b), we decided to implement the code in the Julia programming language due to its efficiency and high performance [12]. Figure 3 displays the computed amplitudes for the classical and relativistic case over the time interval *<sup>t</sup>* <sup>∈</sup> [0, 50] in SI units of seconds, again with all parameters set equal to unity, except for the spring constant *<sup>k</sup>* = <sup>2</sup> kg/s<sup>2</sup> . Using the Julia library

DifferentialEquations.jl [13], we chose the algorithm based on first-order interpolation with at least stepsize <sup>Δ</sup>*<sup>t</sup>* = 0.05 <sup>s</sup> and a relative tolerance of 10−<sup>8</sup> in the adaptive timestepping. Note that a long-time analysis has shown that no significant numerical errors materialize.

As shown in Figure 3, the difference between the simple classical result, resulting from Equation (3) with conditions Equation (8),

$$\mathbf{x}(t) = \frac{m\_0 \mathbf{g}}{k} \left[ \cos \left( \sqrt{\frac{k}{m\_0}} t \right) - 1 \right] \tag{25}$$

and the relativistic numerical estimate is quite pronounced—not only in amplitude but also in phase, with a positive phase shift. The physical results for the relativistic oscillator always possess longer amplitude and period than the classical analogue by cause of the relativistic mass increase.

**Figure 3.** Amplitudes for the relativistic harmonic oscillator compared to the classical model, both suspended in a uniform gravitational field, either using the simple classical solution, viz. Equation (3), or the relativistic estimates resulting by numerical integration of Equations (19a) and (19b). The gravitational field is strong with *<sup>g</sup>*/*c*<sup>2</sup> <sup>∼</sup> <sup>1</sup> <sup>m</sup>−1, choosing again all parameters unity, except for *<sup>k</sup>* = 2 kg/s2, see Figure 2. The two initial conditions are *<sup>x</sup>*(0) = *<sup>x</sup>*˙(0) = 0, see Equation (8).

For an extended time interval *<sup>t</sup>* <sup>∈</sup> [0, 100] in units of seconds, Figure <sup>4</sup> shows the difference between relativistic and classical predictions, <sup>Δ</sup>*x*(*t*) = *xrel*(*t*) <sup>−</sup> *xclas*(*t*), as already individually shown in Figure 3. As found before from the two curves in Figure 3, the relativistic corrections—now more easily recognized—propagate with a positive phase while significantly modulating the amplitude of the classical model. These contributions are substantial and cannot be neglected. Remarkably, these relativistic corrections take the shape of pronounced wave packets.

Figure 5 depicts the phase space for the relativistic harmonic oscillator with potential *Vg* in comparison with the classical model, using the same parameters as before, viz. Figure 2. This phase portrait provides a global overview about the dynamics of the oscillating system and shows the quantitative difference between the two models.

The relativistic phase-space trajectories emerge from the equation of motion, Equation (18), with the preserved quantity

$$\frac{m\_0 c^2}{\sqrt{1 - \frac{\chi^2}{c^2}}} + \frac{1}{2} k x^2 + V\_\mathcal{g}(\mathbf{x}) = E\_0 + m\_0 c^2 \tag{26}$$

which represents energy conservation with total energy *<sup>E</sup>*<sup>0</sup> <sup>+</sup> *<sup>m</sup>*0*c*<sup>2</sup> in Equation (7). In the non-relativistic limit, we may use *<sup>γ</sup>* <sup>≈</sup> <sup>1</sup> + <sup>1</sup> 2 *x*˙ 2 *<sup>c</sup>*<sup>2</sup> and *Vg*(*x*) <sup>≈</sup> *<sup>m</sup>*0*gx* in Equation (26) to obtain

$$
\frac{1}{2}m\_0\dot{\mathbf{x}}^2 + \frac{1}{2}k\mathbf{x}^2 + m\_0\mathbf{g}\,\mathbf{x} = E\_{0\prime} \tag{27}
$$

which is just the well-known result for the harmonic oscillator with uniform Newtonian gravity.

The classical solutions are then reproduced in the phase plane (*x*, *<sup>x</sup>*˙) as the level curves of Equation (27) by varying *E*<sup>0</sup> ≥ 0. Similarly, we obtain the relativistic solutions in the phase plane by drawing the level curves of Equation (26). With the chosen initial conditions in Equation (8), we have *<sup>E</sup>*<sup>0</sup> <sup>=</sup> 0 for both cases (which is equivalent to *<sup>E</sup>* <sup>=</sup> *<sup>m</sup>*0*c*2). Figure 5 shows the two corresponding curves. Observe that Equations (26) and (27) always produce phase-space trajectories (relativistically and classically) which are closed. As a consequence, both types of solutions have to be periodic. The phase portraits of both cases represent center stable dynamics. Another characteristic effect is that the phase-space path for the relativistic case is larger than for the classical path—an effect which also has been observed for the relativistic pendulum [2] and a harmonic oscillator satisfying a relativistic isochronicity principle [4].

**Figure 4.** Amplitude changes, <sup>Δ</sup>*x*(*t*) = *xrel*(*t*) <sup>−</sup> *xclas*(*t*), characterizing the relativistic corrections for the classical harmonic oscillator in a strong uniform gravitational field, corresponding to the physical configuration of Figure 3, but for the larger time interval *<sup>t</sup>* <sup>∈</sup> [0, 100] in units of seconds. These corrections are significant and modulate in both amplitude and phase.

**Figure 5.** The phase-space diagrams for the harmonic oscillator in a uniform gravitational field corresponding to the classical and the relativistic case. We are using the same parameters as in Figure 2. The characteristics for both cases are similar: the closed trajectories represent periodic solutions, and their phase portraits represent center stable dynamics.

#### **5. Conclusions**

We have derived and studied the relativistic generalization of an oscillating body in motion under a Hookean potential and in parallel alignment with a uniform gravitational field.

If gravity is strong, these relativistic corrections differ substantially from the classical predictions and are relevant. These amplitude corrections can be pictured as fluctuating wave packets traveling on top of the classically predicted oscillations.

Quantitatively, the relativistic oscillator always acquires longer amplitudes and periods than the classical analogue. This is the dynamical effect due to the increase in mass of a moving object as dictated by special relativity. In spite of this fundamental difference between the relativistic and non-relativistic framework, we have proven that the corresponding relativistic equation of motion still maintains unique periodic solutions, similar to the well-known classical case.

For practical purposes, we have also presented approximations in the non-relativistic limit (keeping the next-to-leading-order term) for the relativistic gravitational potential and for the equation of motion of the harmonic oscillator—an approximate equation which could still be solved analytically in terms of elliptic integrals. However, all estimates in the fully relativistic regime have to be solved by numerical integration. Toward this end, we implemented the integration code with high precision in the Julia programming language by using efficient first-order interpolation. This simulation also allowed to represent the classical and relativistic phase space and to explore the dynamics of both models with their common and quantitatively different features.

**Funding:** This work has been supported by the Spanish *Ministerio de Economía y Competitividad*, the European Regional Development Fund (ERDF) under grant TIN2017-89314-P, and the *Programa de Apoyo a la Investigación y Desarrollo 2018* (PAID-06-18) of the Universitat Politècnica de València under grant SP20180016.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflicts of interest.

### **References**

