## Abstract

In this paper we report on the feasibility of light confinement in orbital geodesics on
stationary, planar, and centro-symmetric refractive index mappings. Constrained to fabrication and
[meta]material limitations, the refractive index, *n*, has been
bounded to the range: 0.8 ≤ *n*(*r⃗*) ≤ 3.5.
Mappings are obtained through the inverse problem to the light geodesics equations, considering
trappings by generalized orbit conditions defined *a priori*. Our simulation results
show that the above mentioned refractive index distributions trap light in an open orbit manifold,
both perennial and temporal, in regards to initial conditions. Moreover, due to their
characteristics, these mappings could be advantageous to optical computing and telecommunications,
for example, providing an on-demand time delay or optical memories. Furthermore, beyond their
practical applications to photonics, these mappings set forth an attractive realm to construct a
panoply of celestial mechanics analogies and experiments in the laboratory.

© 2013 OSA

## 1. Introduction

Einstein’s general relativity theory allows us to understand the dynamics of mass and
massless particles as they travel through space-time. From it, we learn that matter and geometry are
inextricably related; the latter transfigures the metric of space-time, and this alteration in turn
modifies the geodesics, the paths of particles, resulting in a manifold of relativistic effects,
*v. gr.* the apparent shift in the position of the stars due to the
*gravitational lensing effect*, depicted in Fig.
1(a), where photons traveling near a massive star are deflected from their otherwise straight
path due to the curvature of space-time; or trapping in multifarious orbits around black holes with
respect to the initial conditions and *shape* of the gravitational attractor
[1–3], see Figs. 1(b) and 1(c).

Keeping in mind the previous illustrations, let us momentarily examine a more mundane experiment. Imagine a photon moving in a region of space filled with some monotonically increasing refractive index. What we observe as the photon moves in the medium is that its path describes a curve like the one depicted in Fig. 1(d). If the gradient of the refractive index is tailored properly, it would be possible to build an optical analogy to gravitational bending (see Fig. 1(a)) for example. This leads us to wonder whether it is possible to fabricate an optical analogue to other sidereal geodesics by means of a differential control on the index of refraction.

The latter analogy reveals a connection between both examples, a weft where both realms
[sidereal and laboratory] can be understood as manifestations of the same
mathematical ground, albeit significantly different physics, *differential geometry*.
Einstein’s general theory of relativity tells us that geometry and gravity are linked by
means of the metric tensor of space-time; in optics the intricacy of the relation between geometry
and the physical properties of materials (permittivity and permeability) has been evinced by
miscellaneous methods, most based on differential geometry [4–7].

These methods can be used to show that the optical example of differential light bending can be understood in terms of geometry, essentially by bending the space on which the light propagates, hence the metric is related to the electromagnetic tensors [4, 6]. This realization paved the way for the theoretical realization of optical cloaks, super-lenses, optical illusion devices and recently black-hole analogues in moving refractive index perturbations [6–14].

Concerning trapping and light confinement in particular, different methods have been proposed [2, 15–17]. Of particular interest is Genov et al study, which presented an isotropic inhomogeneous refractive index mapping which was able to confine light mimicking a closed orbit around a black hole [15]. As the authors noted, such mapping could be implemented on a dielectric photonic crystal, composed of air-gaps in Gallium-Arsenide.

In our work, we posit theoretical traps that provide an optical analogue to celestial mechanics to realize light confinement, trapping in a region of space, by means of two-dimensional stationary refractive index mappings, which could be implemented under current technological and photo-refractive [meta]materials, constraint to optical frequencies, which could be easily integrable on all-optical-chips [18–21]. The mathematical and physical background to make these effects possible could bring forth an exciting ground to test sidereal mechanics in the laboratory; and provide the key to enable a miscellany planar optical systems that are of great interest for diverse photonic applications, ranging from time delays, to temporal optical memories, or random resonators, to name a few; thence its significance.

## 2. Deriving the model

We regard a model equation based in the paraxial approximation of non-linear electrodynamics, by
looking into a medium where the index of refraction can be modulated in space and time by some
smooth varying function *f*, e.g. in photo-refractive media index of refraction is a
real valued function *f* of the incident light intensity *I*;
therefore *n*(*x⃗*) =
*f*(*I*(*x⃗*)), where
*x⃗* is the space-time quadrivector which corresponds to the coordinates
{*x*^{0} : *x*^{1},
*x*^{2}, *x*^{3}}. In this representation
*x*^{0} corresponds to the time coordinate, and the rest to the usual
three-dimensional space. In this context the photon’s path is governed by the geometry of
the space the light is propagating on, which in general has the metric tensor
*g _{ij}*, and the infinitesimal length element

*ds*

^{2}=

*g*. Therefore, if we parameterized the path followed by light with reference to the proper time

_{ij}dx^{i}dx^{j}*τ*, i.e.

*x⃗*(

*τ*), we can analytically derive the path equations by means of Lagrange-Euler equations, whit a Lagrangian of the form:

*τ*, i.e. the

*proper time*.

The physical description behind Eq. (1) is
relevant to optics, as we set forth earlier in the introduction. If the metric were indeed related
to the electromagnetic tensors, permittivity and permeability, then it could be possible to write
the light-path in the paraxial approximation as a function of either the refractive index $n=f\left({x}^{0},\overrightarrow{r}\right)=\sqrt{\epsilon \left(\overrightarrow{r},t\right)\mu \left(\overrightarrow{r},t\right)}$, or the metric *ds*. Hence, it could be possible to use
Fermat’s principle, where Eq. (1) would take
the role of an *optical Lagrangian*, i.e. the solution to the variational problem
which states that given a refractive index *n*(*x*^{0},
*r⃗*) the path light follows is given by the extrema of:

*j*= 0,..., 3, and

*δ*indicates the functional maximization. Recursively, the solution to Eq. (2) is the Lagrange-Euler Eq. (1). Meanwhile, observe that in the integral of Eq. (2) the term

*ds*takes into account the metric by virtue of the covariant metric tensor

*g*of the curved space:

_{ij}Yet, we are still to determine the relation between the metric tensor
*g _{ij}* and the refractive index

*n*. This can be achieved by means of conformal mappings, where it is possible to build the relation between the metric and the electromagnetic properties of a [meta]material [4–6]. Specifically, if the line element of space-time takes the form

*ds*

^{2}=

*g*

_{00}

*dx*

^{0}

*dx*

^{0}−

*g*, it is possible to map the curved space-time to an isotropic, non-dispersive and non-absorbing medium as [4, 15]:

_{ii}dx^{i}dx^{i}*g*= det(

*g*). Furthermore if we consider Minkowski’s metric of special relativity, the tensor

_{ij}*g*has the entries:

*n*=

*n*(

*x*

^{0},

*r⃗*) can be found using Eq. (4). And consequently, the length element

*ds*of Eq. (2) can be written as:

#### 2.1. Path equations and the inverse problem

Accordingly, from the Lagrange-Euler (1) we can derive the light-path, which in cylindrical coordinates takes the form:

*τ*, e.g.

*ṙ*=

*∂r*/

*∂τ*. Equation (7) through Eq. (9) describe the parametrized path of a light beam traveling in a space-time with such geometry or equivalent refractive index.

We will describe light being trapped, when the geodesics are perpetually bounded to some region of space, for as long as the material properties remain unchanged. In this sense a light trapping geodesic is one where any of the following conditions arise:

*ṙ*,*ϕ̇*→ 0,*ṙ*→ 0 and*ϕ*≤_{min}*ϕ*≤*ϕ*,_{max}*ϕ̇*→ 0 and*r*≤_{min}*r*≤*r*,_{max}*r*≤_{min}*r*≤*r*and_{max}*ϕ̇*≤_{min}*ϕ̇*≤*ϕ̇*,_{max}

*τ*≥

*τ*

^{*}; where

*τ*

^{*}marks the moment at which the initial conditions for trapping are met. The first condition is only possible if the refractive index

*n*goes to infinity, or if the refractive index co-moves with the wave [16, 17]. The other three cases can be realized in

*stationary*refractive index [15, 22]. Thus, they can be imposed as general conditions to solve the inverse problem to Eq. (7) through Eq. (10).

The inverse problem begins by defining a set of *a priori* conditions that
describes a generalized trapping scenario. If we impose premises 1 to 4, described earlier, the path
equations reduce to new set of equations where the affine parameter is time itself; whence we can
rework Eq. (7) through Eq. (10) into polar coordinates:

Since we are interested in the general stationary case, we impose the fourth condition above,
ensuing 0 ≤ *ṙ* ≤ *ṙ _{max}*
and

*ϕ̇*≤

_{min}*ϕ̇*≤

*ϕ̇*, on Eq. (11) and Eq. (12), and find a non-singular refractive index distribution,

_{max}*n*(

*r*,

*ϕ*), that satisfies them. Notice that we require the refractive index to be smoothly varying and

*C*

^{2}continuous, i.e. ∇

*n ·e⃗*≤

^{i}*M*, where

*M*is a real valued number, and

*e⃗*a vector in a give direction in space-time. The general solution to this problem, however, is complex and a unique solution is not guaranteed [18, 23, 24].

^{i}The latter being said, the problem can be simplified further if we recall that a trapping orbit
occurs when the radius *r* is bounded to a region of space, i.e.
*r _{min}* ≤

*r*≤

*r*, for all

_{max}*ϕ*; whence

*ṙ*(

*r*) =

_{max}*ṫ*(

*r*) = 0, and

_{min}*r*can be written in terms of the angle

*ϕ*or the time

*t*. Hence a natural solution to the inverse problem is an autonomous refractive index that varies with the radial variable, i.e.

*n*=

*n*(

*r*). This refractive index map resembles the centro-symmetric characteristics of gravitational fields, where, as pointed out independently by Genov and Ni [15, 22], light would behave similarly to that traveling in curved space-time. Such distribution can be constrained to the conditions listed previously.

This result is significant in terms of fabrication where novel photo-refractive [metal]materials could be used. The refractive index in these composites can be modulated by the intensity of the change-inducing beam; which can be controlled by shaping its transversal profile [18–21]. Accordingly several families of refractive index distributions can be achieved, e.g. Gaussian, Laguerre, or Bessel, which are all centro-symmetrical.

#### 2.2. System stability and dynamic equillibrium

Since the refractive index is a function of the radius, Eq. (11) and Eq. (12) can be rewritten in the form (13):

*ℓ*is a constant of integration. This equation is reminiscent of orbital mechanics, with the sole exception that

*μ*, the mass in celestial mechanics, varies with the radius. Substituting the latter result into Eq. (14), and rearranging yields:

*ratio of external energy supply*, which is a measure of the equilibrium stability of the system [25, 26]; simplifying we get:

By virtue of Bertrand’s theorem [15, 27], we know that *ṙ* ≠ 0,
for all time *t*; then from premises 3 and 4 it follows that
*ṙ _{min}* ≤

*ṙ*≤

*ṙ*, where

_{max}*ṙ*< 0; and from the definition of

_{min}*n*(

*r*) as a monotonous decreasing function it follows that

*∂*ln

_{r}*n*≤ 0. Therefore, as the radius decreases, the radial velocity goes to its minimum and

*ℰ*≤ 0, therefore the system approaches the equilibrium point and the amplitude of the path decreases. On the other hand, as the radius increases, the velocity tends to a maximum, and

_{t}*ℰ*≥ 0, thus the system is driven away from the equilibrium point and the amplitude of the path increases. In Fig. 2 we show

_{t}*ℰ*for a Gaussian like refractive index described in section 3.1. Observe that in the light confinement scenario we will see an oscillatory behavior in

_{t}*ℰ*.

_{t}To conclude we study the stability of the orbits. by means of Lyapunov’s stability theory for the governing Eq. (17), and derive the constraints applicable to the refractive index mapping to attain stable orbital paths [15]. First, let us rewrite Eq. (17) as:

*β*=

*ϕ̇*and

*χ*(

*r*) =

*∂*ln

_{r}*n*. Then, if we introduce small perturbations

*r*=

*r*+

_{o}*δ*and ϰ = ϰ

*+*

_{o}*ε*, expand

*χ*(

*r*) in a Taylor series around

*r*, and drop all higher order terms for

_{o}*δ*and

*ε*; we can rewrite Eq. (19) as a linear system of equations for the perturbation terms:

*χ*and

_{o}*χ*′

*are the linear terms of the Taylor expansion of*

_{o}*χ*(

*r*) around

*r*. Ensuing the eigenvalues of the matrix above are given by:

_{o}The geodesics will be stable if the real part of the *eigenvalues* of Eq. (20),
ℜ{*λ*_{1,2}}, are negative or equal to zero.
Recall that 0 ≤ *β* ≤
lim_{r}_{→∞}*n*(*r*)^{−1},
and *χ*(*r*) ≤ 0; and therefore the condition for
Lyapunov stability is reached if:

This constrained is met for all points *r* near the stability point
*r _{o}* = 0, provided the radial velocity meets

*ṙ*≤ 0. In general, since the refractive index is monotonous decrescent, the system will be Lyapunov stable at points where

_{o}*ṙ*≤ 0, see, for example, the phase-space plots in Figs. 6(b) and 8(b). In addition, it is interesting to analyze whether the system will be Lyapunov stable under the conditions necessary to sustain a circular orbit. In this case

*ϰ*= 0, and

_{o}*r*=

_{o}*a*, constant, for all time

*t*; consequently

*V*= 0,

*W*≤ 0, and the second term in the equality (22) is imaginary, provided that

*χ*(

*r*) ≤ 0; thus resulting in ℜ{

*λ*

_{1,2}} =

*V*= 0. These conditions, however, result in a different family of refractive indexes. In our study we focus on Gaussian-like refractive index distributions, feasible through photo-refraction, where

*ϰ*≠ 0, ergo the system will not describe closed circular orbits.

In summary, light confinement will occur as long as the conditions cited earlier apply. Under
them the system will oscillate around the equilibrium point as depicted by the *ratio of
external energy supply*. In the next section we will present the confinement geodesics
resulting from these refractive index distributions. We will show that the geodesics in this system
mimic those of a spinning black hole and other celestial objects, thus reproducing the paths of
light traveling in curved space.

## 3. Results and discussion

In this section we explore three refractive index distributions based on Gaussian profiles. The first two mimic light trapped in the vicinity of a black hole, and the third case that of a particle moving in-between two planets.

#### 3.1. Guassian distribution, single attractor

A Gaussian refractive index distribution has the form:

Then making it possible to constrain the maximum and minimum of the refractive index distribution
to feasible values, e.g. 0.8*n* ≤ 3.8. Hence, given a range of values for the
back-ground and maximum refractive index, *n _{c}* and

*n*, respectively, we can evaluate different light confinement scenarios, see Fig. 3.

_{a}As discussed in the previous section the orbital paths are open, while their shape and stability
depend significantly on the initial conditions, *n*(*r*) and
{*r _{o}*,

*ṙ*,

_{o}*ϕ*,

_{o}*ϕ̇*}. As portrayed by Fig. 3 different initial conditions will result in new orbit manifolds with distinct characteristics, albeit maintaining a common morphology. As it can be seen from these figures the orbital paths are highly sensitive to the modifications in the refractive index and the initial conditions; a change in ∼ 0.5% in

_{o}*n*alters the orbit, as seen from the plot of radii extrema in Fig. 4.

_{c}In Fig. 5 we show the 2D open orbit for a refractive index
map where *n _{c}* = 0.8 and

*n*= 3.0. As detailed in Fig. 6(a), the radial velocity describes an oscillatory form, which indicates that the radius follows a similar motion, swinging between two extrema, as expected in a confine and open orbital motion. The latter can be observed in greater detail by analyzing the phase-space, plotted in Fig. 6(b); notice that when the velocity

_{a}*ṙ*reaches zero the orbit radius

*r*reaches an extrema; for a closed single loop the amplitude of the velocity oscillation would be zero, and the phase space would become a horizontal line at

*ṙ*= 0. The behaviour of the angle versus the radial velocity, and consequently to the radius movement, is described by Fig. 6(c), where the absolute value of the radial velocity is plotted against the angular position, which results in a periodic motion. Observe that in the limit at every angle there would be a time at which the radial velocity would be zero, which is typical of open orbits.

#### 3.2. Mexican hat distribution

We study a variation on the Gaussian trap to understand the potential modifications on the
trapping orbits when large modifications ∼ 20% of the refractive index take place.
The mathematical description of the *mexican hat* is:

Using the same conditions as in the previous analysis we set *n _{c}*
= 0.8,

*n*= 3.0, and set

_{a}*n*= 0.2. The results, shown in Figs. 7 and 8, display an orbit whose appearance resembles that of the regular Gaussian trap studied earlier. However, a closer look reveals that the trap confinement area is ∼ 2% slimmer, and the oscillatory frequency of the radial and angular velocities, Fig. 8(a), increases by ∼ 30%. Whereas the radial phase-space, Fig. 8(b), shows a smaller eccentricity compared to its Gaussian counterpart, and the absolute value of the radial velocity versus the angle, Fig. 8(c), describes the typical behaviour of an open orbit.

_{d}These variations in the trapping conditions could potentially be used to differentiate trapping in transient optical memories to allocated binary data.

#### 3.3. Double attractors

Interested in observing phenomenologically trapping in a binary system we evinced a planar refractive index map resulting from the superposition of two Gaussian functions,

The resulting orbits are shown in Fig. 9. A binary system,
however, presents a relatively higher sensitivity to modifications in the refractive index, where a
distortion of Δ*n* > 0.03 modifies the trapping orbits, in some cases
releasing the light beam.Yet, it is possible to confine the light beam to an orbit which mimics the
behaviour of Newton’s three body problem, where, as expected, small perturbations result in
sundry possible solutions to the path equation, i.e. chaos.

#### 3.4. The space-time equations, trapping

Briefly, we show that the same Gaussian solution explored in section 3.1 fulfills the trapping conditions in the full eikonal space-time, Eq. (7) through Eq. (9). The results are shown in Fig. 10.

In this case we observe a similar behaviour to that of the Gaussian trap in the stationary approximation, the difference in the amplitude values is due to normalization. An oscillatory behaviour in the radius and angular velocities is described in Figs. 10(b) and 10(c); as it would be expected from a centro-symmetrical refractive index mapping.

#### 3.5. Material realization

Recent advancement in photo-refractive materials [18–21] could enable the construction
of the light confinement devices discussed. For example, a prototype Gaussian trap could be
fabricated by stacking different photo-refractive [meta]materials like
*As _{x}S*

_{1−x}in multiple layers. Recalling that in photo-refraction the optical properties of the material, permeability and permittivity, change as a function of the light intensity, allowing the fabrication of a continuously varying refractive index as the one required by a centro-symmetrical attractors. Moreover, since photo-refractivity is elastic,

*i.e.*reversible, it allows one to modify the refractive index of the layers combined in time, thence enabling or disabling the trap at will.

## 4. Conclusions

We have demonstrated that, under the paraxial description of light in heterogeneous media, it is possible to use centro-symmetrical refractive index distributions as attractors for light, which are able to perform light trapping or confinement in open orbits. These mappings are bound within fabrication constrains, and due to the paraxial approximation they are several times larger than the mean beam width and wavelength, ∼ 50 times larger. The proposed mappings have potential applications to transient optical memories, delays, concentrators, random cavities, and beam stirrers, to cite a few; which could enable the next generation of photonic integrated circuits, PICs, processors and computation systems. Due to their sensitivity to the change in the refractive index these devices could be employed as a sensor platform. Furthermore, since they are drawn in close analogy to celestial phenomena they could be used as laboratory setups to test a large variety of effects in celestial mechanics.

## Acknowledgments

The authors wish to acknowledge David Ketcheson, Carlos Argaez, Pedro Guerra, and Jose Alcantara-Felix for their suggestions, corrections and valuable discussions while preparing this manuscript.

## References and links

**1. **J. Levin and G. Perez-Giz, “A periodic table for black hole
orbits,” Phys. Rev. D **77**, 103005 (2008) [CrossRef] .

**2. **F. T. Hioe and D. Kuebel, “Characterizing planetary orbits and trajectories of light
in the Schwarzchild metric,” Phys. Rev. D **81**, 084017 (2010) [CrossRef] .

**3. **T. Müller and J. Frauendiener, “Studying null and time-like geodesics in the
classroom,” Eur. J. Phys. **32**, 747–759 (2011) [CrossRef] .

**4. **U. Leonhardt and T. G. Philbin, “Transformation optics and the geometry of
light,” in *Progress in Optics*53, 69–153, E. Wolf, ed. (Elsevier, 2008) [CrossRef] .

**5. **H. Chen, C. T. Chan, and P. Sheng, “Transformation optics and
metamaterials,” Nat. Mater. **9**, 387–396 (2010) [CrossRef] [PubMed] .

**6. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic
fields,” Science **312**, 1780–1782, 2006 [CrossRef] [PubMed] .

**7. **D. Schurig, J. B. Pendry, and D. R. Smith, “Calculation of material properties and ray tracing in
transformation media,” Opt. Express **14**, 9794–9894 (2006) [CrossRef] [PubMed] .

**8. **U. Leonhardt and Piwnicki, “Optics of non-uniformly moving
media,” Phys. Rev. A **60**, 4301–4312 (1999) [CrossRef] .

**9. **U. Leonhardt and Piwnicki, “Relativistic effects of light in moving media with
extremely low group velocity,” Phys. Rev. Lett. **84**, 822–825 (2000) [CrossRef] [PubMed] .

**10. **I. Brevik and G. Halnes, “Light rays at optical black holes in moving
media,” Phys. Rev. D **65**, 024005 (2001) [CrossRef] .

**11. **R. Liu, C. Ji, J. J. Mock, J. Y. Chin, T. J. Cui, and D. R. Smith, “Broadband ground-plane cloak,”
Science **323**, 366–369 (2009) [CrossRef] [PubMed] .

**12. **T. Ergin, N. Stenger, Brenner, J. B. Pendry, and M. Wegener, “Three-dimensional invisibility cloak at optical
wavelengths,” Science **328**, 337–339 (2010) [CrossRef] [PubMed] .

**13. **Y. Lai, J. Ng, H. Chen, D. Han, J. Xiao, Z.-Q. Zhang, and C. Chan, “Illusion optics: the optical transformation of an object
into another object,” Phys. Rev. Lett. **102**, 253902 (2009) [CrossRef] [PubMed] .

**14. **D.-H. Kwon and D. H. Werner, “Flat focusing lens designs having minimized reflection
based on coordinate transformation techniques,” Opt.
Express **17**, 7807–7817 (2009) [CrossRef] [PubMed] .

**15. **D. Genov, S. Zhang, and X. Zhang, “Mimicking celestial mechanics in
metamaterials,” Nat. Phys. **5**, 687–692 (2009) [CrossRef] .

**16. **Piwnicki and U. Leonhardt, “Optics of moving media,”
App. Phys. B **72**, 51–59 (2001) [CrossRef] .

**17. **S. L. Cacciatori, F. Belgiorno, V. Gorini, G. Ortenzi, L. Rizzi, V. G. Sala, and D. Faccio, “Spacetime geometries and light trapping in travelling
refractive index perturbations,” New J. Phys. **12**, 095021 (2010) [CrossRef] .

**18. **D. P. San-Román-Alerigi, T. K. Ng, M. Alsunaidi, Y. Zhang, and B. S. Ooi, “Generation of *J*_{0}-Bessel-Gauss
beam by an heterogeneous refractive index map,” J. Opt. Soc.
A **29**, 1252–1258 (2012) [CrossRef] .

**19. **A. Zakery, “Optical properties and applications of chalcogenide
glasses: a review,” J. Non-Crystaline Solids **330**, 1–12 (2003) [CrossRef] .

**20. **S. B. Kang, M. H. Kwak, B. J. Park, S. K., H. C. Ryu, D. C. Chung, S. Y. Jeong, D. W. Kang, S. K. Choi, M. C. Paek, E. J. Cha, and K. Y. Kang, “Optical and dielectric properties of chalcogenide glasses
at terahertz frequencies,” ETRI J. **31**, 667–674 (2009) [CrossRef] .

**21. **D. P. San-Román-Alerigi, D. H. Anjum, Y. Zhang, X. Yang, A. B. Slimane, T. K. Ng, M. N. Hedhili, M. Alsunaidi, and B. S. Ooi, “Electron irradiation induced reduction of the permittivity
of chalcogenide glass (As_{2}S_{3}) thin fim,” J.
App. Phys. **113**, 044116 (2013) [CrossRef] .

**22. **X. Ni and Y.-C. Lai, “Transient chaos in optical
metamaterials,” Chaos **21**, 033116 (2011) [CrossRef] [PubMed] .

**23. **B. Gangand and L. Peijun, “Inverse medium scattering problems for electromagnetic
waves,” J. App. Math. **65**, 2049–2066
(2005).

**24. **M. Piana, “On uniqueness for an isotropic inhomogeneous inverse
scattering problems,” Inv. Problems **14**, 1565–1579 (1998) [CrossRef] .

**25. **D. W. Jordan and P. Smith, *Nonlinear ordinary differential equations* (Oxford
University Press, 2007).

**26. **F. Verhulst, *Nonlinear differential equations and dynamical systems*
(Springer Verlag, 1990) [CrossRef] .

**27. **Y. Zarmi, “The Bertrand theorem revisited,”
Am. J. Phys **70**, 446–449 (2002) [CrossRef] .