Single Particle Closed Orbits in Yukawa Potential

Rupak Mukherjee†123 , Sobhan Sounda‡3

arXiv:1705.02444v1 [[Link]-ph] 6 May 2017

1. Institute for Plasma Research, HBNI, Gandhinagar, Gujarat, India.

2. Ramakrishna Mission Vivekananda University, Belur Math, Howrah, Kolkata, India.

3. Ramakrishna Mission Residential College (Autonomous), Narendrapur, Kolkata, In-

dia.

† rupakmukherjee01@[Link]
‡ sounda6@[Link]

May 9, 2017

Abstract:

We study the orbit of a single particle moving under the Yukawa potential and observe
the precessing ellipse type orbits. The amount of precession can be tuned through the cou-
pling parameter α. With a suitable choice of the coupling parameter; we can get a closed
bound orbit. In some cases we have observed some petals which can also have a closed
bound nature with an appropriate choice of the coupling constant. A threshold energy has
also been calculated for the boundness of the orbits.

Keywords:
Yukawa Potential, Precessing Ellipse, Closed-Bound Orbit, Critical Point

PACS No:
45.20.D, [Link], [Link], [Link]
1 Introduction:
It is widely known that there are only two types of central potential 1r and r2 in which all
finite motions take place in closed paths[1]. Some exceptions of the statement above have
been reported recently[2, 3].
On the other hand Yukawa potential or screened coulomb potential of the form V (r) =
r
− αr e− λ , ( λ is the screening parameter that determines the range of this interaction)
has a long standing legacy to represent various physical systems. Its application ranges
from astronomy, high-energy physics, nuclear physics, condensed matter physics to plasma
physics and many other branches of physics. In the domain of plasma physics this potential
is used mostly in Strongly Coupled Plasmas which has several applications in modelling
the cores of white dwarfs, planetary rings, atmospheric lightning, molten salts and plasma
technology. Hence quite naturally one may ask what is the orbit of a single particle under
this screened coulomb potential. Can we expect some bound orbits for the non-zero value
of α?
This study is a more general case of the previous studies [2, 3] that has been done till now
regarding the closed bound orbits distinct from coulomb or harmonic oscillator potentials
because by the truncation of the exponential series in the Yukawa potential one gets back
the potentials assumed in the literature so far. Hence this study is a superset of the previous
studies performed so far to the best of our knowledge.
The paper is organised as follows. Section 2 discusses the theory of central force motion for
any arbitrary potential. Following the theory we have calculated the trajectories for Yukawa
potential in section 3 and the thrust has been given on the calculation for threshold energy
for bound orbit in Yukawa potential in section 4. In section 5, some possible applications
of this study has been mentioned.

2 A Brief Discussion of Central Force Motion

The total energy of a particle moving under a central force is given by,

1 2 J2
mṙ + + V (r) = E, J = mr2 θ̇
2 2mr2
where m is the mass of the particle, r is the radial distance of the particle from the force
center, J is the angular momentum of the particle, V is the potential of the particle and θ
is the angular coordinate of the particle with respect to some reference axis. For r = u1 i.e.
J du
ṙ = − m dθ
; the above equation looks like
 2
1 J 2 du J 2 u2
m 2 + + V (u) = E
2 m dθ 2m
dE
For energy to be constant, i.e. = 0,
dθ

d2 u J 2 u2
 
m d m dVef f
2
=− 2 V + =− 2
dθ J du 2m J du

2
where dJdθ
= 0 since, for central force the angular period ∆θ of the radial oscillators r(θ)
is independent of the angular momentum J. In dimensionless coordinates, u = ua and
V = Jm2 a2 V the equation of motion becomes,

d2 u dVef f
2
=−
dθ du
with Vef f = V (u) + 12 u2 . Now if the form of the potential is known, we can find the
trajectory of the particle. The threshold energy value for a bound orbit in any central
force can be calculated from the total energy (h) conservation relation,
Zu
1 2 du
my + Vef f (u) = h, where, Vef f (u) = − F (u)du and y =
2 dθ
0
r
2
⇒y=± [h − Vef f (u)] (1)
m
Now, if we differentiate the above equation we get,
dy 1 dVef f
y = (2)
du m du
dV ef f dy
Here, for du = 0, if du is not undefined, i.e. the curve plotted from equation (1), does
not cut the u axis perpendicularly, then y = 0 which in turn represents h = Vef f i.e. the
minimum threshold value of energy for bound orbit.
dy
Now if du is undefined, i.e. the curve plotted from equation 1, cuts the u axis perpendic-
dVef f
ularly, then from the roots (uc , 0) or the critical points of the equation du = 0 we find
dV (uc ) dVef f
some equilibrium points of the system. Now from the relations du = 0 and du =0
we can conclude that V (u) has either a relative extremum or a horizontal inflection point
at u = uc .
If,
I. the Vef f (u) has a relative minimum at u = uc , then Vef f (uc ) = h0 (threshold or minimum
value for energy) and the critical point is a centre and is stable.
II. the Vef f (u) has a relative maximum at u = uc , then the critical point is a saddle point
and is unstable.
III. the Vef f (u) has a horizontal inflection point at u = uc , then the critical point is of a
degenerate type called a cusp and is unstable.

3 Motion of a Single Particle under Yukawa Potential

The Yukawa or screened coulomb potential is given by:
α r
V (r) = − e− λ (3)
r
1 mα 1
⇒ Vef f (u) = u2 − 2 ue− u (4)
2 J a

3
where, λ is the range of the Yukawa potential.
The equation of motion is given by

d2 u
 
1
−u 1 mα
= −u + αe 1+ , where, α = (5)
dθ2 u J 2a

which can be rewritten in terms of energy as:

1 2 1 2 1
my + u − αue− u = h
2 q2
1
⇒ y=± 2h − u2 + 2αue− u (6)

where we have chosen m = 1. The exact values of α differ for different physical cases
depending on the nature of application[11, 12, 13, 16].

The Equation (5) is difficult to solve for an exact analytical solution. For planetary
type motion the equation of motion with our parameters becomes,

d2 u
 
1
− λu 1
+ u = 1 + αe 1+ (7)
dθ2 λu

In ref[10] an approximate analytical solution has been obtained by expanding the R.H.S.
of (7) in a Taylor series and truncating it to the second order. Thus the equation looks
like,
 2 
d2 u ap0 −ap0 /λ
 
+u 1−α e = up
dθ2 λ2
ap0 a2p0
  
1 −ap0 /λ
up = 1 + αe 1+ − 2
p0 λ λ

where ap0 ∼ p0 = 1/u0 and u0 is the unperturbed (α = 0) solution. The solution of the
above equation given by,

u(θ) = up + ue cos ω(θ − θ0 )

1/2
ω = 1 − α(a2p0 /λ2 )e−ap0 /λ


is only valid for α << 1 and can deviate significantly for particle orbits which lie away
from the unperturbed (α = 0) value. The aim of our present study is to explore the full
solution space of equation (5) by retaining the complete nonlinear form of the force term.
We obtain these solutions by a numeical solution of the equation.

We adopt another analytical tool viz. a linear stability analysis for a particle under
yukawa type potential. We proceed to some extent and then turn to the numerical analysis
for different α (coupling constant).

4
4 Linear Stability Analysis[7] for Single Particle Mo-
tion in Yukawa Potential
For the sake of analysis, we add a small viscous term proportional to velocity (µ du
dθ
), in the
equation of motion and eventually put µ = 0 for the final calculation.
Thus for Yukawa Potential the modified equation of motion (5) can be written as,
 
00 0 0 1 −1
u (θ) = f (u (θ), u(θ)) = −µu − u + α 1 + e u
u
⇒ u0 (θ) = F (u(θ), y(θ)) = y
 
0 1 −1
y (θ) = G(u(θ), y(θ)) = −µy − u + α 1 + e u
u

Hence, the fixed points f (y ∗ , u∗ ) = 0 are:

(y ∗ , u∗ ) = (0, −0.2114) for α = 0.0005.
(y ∗ , u∗ ) = (0, −0.3046) for α = 0.005.
(y ∗ , u∗ ) = (0, −0.4705) for α = 0.05.
(y ∗ , u∗ ) = (0, −0.7296) for α = 0.5.
(y ∗ , u∗ ) = (0, −0.8094) for α = 1.
The Jacobian for the above set of equation will be given by[4],
 ∂F ∂F   
∗ ∗ ∂u ∂y 0 1
J(u , y ) = ∂G ∂G = α −u 1
∂u ∂y u=u∗ ,y=y ∗ −1 + u3
e −µ u=u∗ ,y=y ∗

The Eigenvalues of the above matrix will be,

√  1 r 
1 2 2
 α −1 
λ1,2 = τ ± τ − 4∆ = −µ ± µ − 4 1 − 3 e u
2 2 u u=u∗

It is checked that the values of ∆ are posi-

tive for all values of α given above. Hence, Value of Delta with alpha
10
for the above parameter values the fixed "Delta" u 1:2

point can not be a saddle point. 8

Now, for ∆ > 0

if τ < 0 and τ 2 − 4∆ > 0 then f (y ∗ , u∗ ) is 6
Delta

a stable node.
if τ < 0 and τ 2 − 4∆ < 0 then f (y ∗ , u∗ ) is 4

a stable spiral. 2

if τ > 0 and τ 2 − 4∆ > 0 then f (y ∗ , u∗ ) is

an unstable node. 0
0.0001 0.001 0.01 0.1 1

if τ > 0 and τ 2 − 4∆ > 0 then f (y ∗ , u∗ ) is alpha

an unstable spiral.
if τ = 0 and τ 2 − 4∆ > 0 then f (y ∗ , u∗ ) is a nutrally stable center.

5
The existance of the limit cycle around each of the fixed points (for different values of α)
has been checked and it is found that for none of the cases there exists any limit cycle.
Hence if we continuously change µ from positive to negative the fixed point changes from
stable to unstable spiral. However at µ = 0 we do not have a true hopf bifurcation because
there are no limit cycles on either side of the bifurcation. This situation is identical to the
case of a damped pendulum or a duffing oscillator.
Further we analyse the case of µ = 0 numerically for the quantitative understanding of the
parameter values with closed orbits..

5 Numerical Analysis for Single Particle Motion in

Yukawa Potential
For numerical analysis, we have used a simple Runge Kutta 4th Order solver and have
used standard accepted algorithm for avoiding the numerical singularities, if encountered.
Interestingly we have shown that for some typical values of α we get a closed bound orbit.
In some cases the orbits have some petals and the number of petals can also be controlled
through the judicious choice of the coupling constant. Here we present some results of our
simulation that helps to understand the dependency on the coupling constant easily.

For α = 0.05; we observe a constant precesion of the ellipse. But the ellipse does not
close itself when it completes one revolution. We have tuned the value of the coupling
constant to 0.0501 and have observed the expected closure.

alpha= 0.05 alpha = 0.0501

1.5 1.5
"alpha_0p05 u 1:2 "alpha_0p0501" u 1:2

1 1

0.5 0.5
u*sin(theta)

u*sin(theta)

0 0

-0.5 -0.5

-1 -1

-1.5 -1.5
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5
u*cos(theta) u*cos(theta)

Figure 1: Precessing Ellipse with α = 0.05 and Closed orbits with α = 0.0501

For α = 0.1; we observe almost same behavior but the amount of precesion has increased.
Still we had to tune the coupling constant to 0.11 to get a closed bound orbit.

6
 alpha= 0.1 alpha = 0.11
1.5 1.5
"alpha_0p1" u 1:2 "alpha_0p11" u 1:2

1 1

0.5 0.5
u*sin(theta)

u*sin(theta)
0 0

-0.5 -0.5

-1 -1

-1.5 -1.5
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5
u*cos(theta) u*cos(theta)

Figure 2: Precessing Ellipse with α = 0.1 and Closed orbits with α = 0.11

For α = 0.5; we observe a small grouping of the precessed orbits. The precession
occurs in a bunch of 2 orbits. Hence we may observe that there are two different values of
precession within one complete revolution. When we set the coupling constant to 0.49 we
suddenly observe a closed bound orbit maintaining the nature of grouping.

alpha = 0.5 alpha = 0.49

1.5 1.5
"alpha_0p5" u 1:2 "alpha_0p49" u 1:2

1 1

0.5 0.5
u*sin(theta)

u*sin(theta)

0 0

-0.5 -0.5

-1 -1

-1.5 -1.5
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5
u*cos(theta) u*cos(theta)

Figure 3: Precessing Ellipse with α = 0.5 and Closed orbits with α = 0.049

For α = 1; we observe two distinct types of orbits, which has been further investigated
in the next step. Here a 3-fold larger orbit as well as a 3-fold smaller orbit are found.
The tuning of the coupling constant has been done accordingly to obtain the closed bound
orbit.

7
 alpha = 1.0 alpha = 1.011
2 1.5
"alpha_1p0" u 1:2 "alpha_1p011" u 1:2

1.5
1

1
0.5

0.5
u*sin(theta)

u*sin(theta)
0

-0.5
-0.5

-1
-1

-1.5
-1.5

-2 -2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5 -1 -0.5 0 0.5 1 1.5 2
u*cos(theta) u*cos(theta)

Figure 4: Precessing Ellipse with α = 1.0 and Closed orbits with α = 1.011

For α = 5; we observe two distinct classes of petals, one within the other. Maintaining
the two distinct petal structures we have been able to find the closed bound orbits for some
parameter value (α).

alpha = 5.0 alpha = 5.1002

10 10
"alpha_5p0" u 1:2 "alpha_5p1002" u 1:2

8 8

6 6

4 4

2 2
u*sin(theta)

u*sin(theta)

0 0

-2 -2

-4 -4

-6 -6

-8 -8

-10 -10
-10 -8 -6 -4 -2 0 2 4 6 8 10 -10 -8 -6 -4 -2 0 2 4 6 8 10
u*cos(theta) u*cos(theta)

Figure 5: Precessing Ellipse with α = 5.0 and Closed orbits with α = 5.1002

For α = 10 also; we observe two distinct classes of petals but the ratio of the bigger to
smaller orbit has changed. Here also we have observed closed orbits.

8
 alpha = 10.0 alpha = 9.91000045528
20 20
"alpha_10p0" u 1:2 "alpha_9p91" u 1:2

15 15

10 10

5 5
u*sin(theta)

u*sin(theta)
0 0

-5 -5

-10 -10

-15 -15

-20 -20
-20 -15 -10 -5 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20
u*cos(theta) u*cos(theta)

Figure 6: Precessing Ellipse with α = 10.0 and Closed orbits with α = 9.91000045528

Thus we observe that increasing the magnitude of coupling constant increases the pre-
cession speed. In some cases it also changes the nature of the orbits. Thus if one can tune
the coupling constant in such a way that after one complete rotation the precessed orbit
matches with the first orbit, it traces back the same path and hence forms a closed bound
orbit. The same is true for the orbits with two distinct types of petals. Thus for Yukawa
potential boundness criteria is just a cautious choice of the coupling parameter that we
have presented extensively for many parameter (coupling constant; α) values in the above
figures.

6 Energy diagrams for yukawa potential:

Now we calculate the threshold energy for the bound orbits. For a constant value of
coupling constant, α = 0.05, we plot equation 6 for different values of energy. (h =
2, 1, 0.5, 0.25, 0.1, 0.05, 0.005). Note that all the curves cut the negative u axis perpendicu-
dVef f
larly representing du to be undefined. This states that we have not reached the threshold
energy value. For energy values less than 0.5 we observe no crossing in the negative u axis
representing that we have gone further below the threshold energy for coupling constant
α = 0.05

9
 alpha = 0.05
2
"h=2" u 1:2
"h=1" u 1:2
"h=0.5" u 1:2
1.5 "h=0.25" u 1:2
"h=0.1" u 1:2
"h=0.05" u 1:2
"h=0.005" u 1:2
1

0.5

0
y

-0.5

-1

-1.5

-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
u

Figure 7: Energy diagram for constant coupling parameter α = 0.05 and varying energy h

Again for a constant h = 1.5 we plot for different values of coupling constant (α =
0.0005, 0.005, 0.05, 0.5, 1)

h = 1.5
2
"alpha=0.0005." u 1:2
"alpha=0.005." u 1:2
"alpha=0.05." u 1:2
1.5 "alpha=0.5." u 1:2
"alpha=1.0." u 1:2

0.5

0
y

-0.5

-1

-1.5

-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3
u

Figure 8: Energy diagram with constant energy h = 1.5 and varying coupling parameter α

Now we calculate the minimum values of energy and the u-nullclines graphically. For
a constant value of α we start with a bound orbit in the negative u axis and keep on
decreasing energy untill the circle reduces to a point. Then if the value of energy is further
decreased; suddenly we observe that the point in the negative u axis disappears. This
corresponds to the minimum value of energy for a closed orbit in the u − θ diagram. In
the figure below we have plotted the minimum value of energy (h) upto which the point in

10
the negative u axis was found. The position of the point gives the value of the u-nullcline
and we have measured that value for different α.

Plot for Minimum Energy

2.5
"alpha=1.0" u 1:2
"alpha=0.5" u 1:2
2 "alpha=0.05" u 1:2
"alpha=0.005" u 1:2
"alpha=0.0005" u 1:2
1.5

0.5

0
y

-0.5

-1

-1.5

-2

-2.5
-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5
u

Figure 9: Threshold energy values for different coupling parameters α

From the above graph we have made a table for coupling constant and minimum thresh-
old energy for a bound orbit (that we have obtained numerically by varying the parameter
h) and the u-nullcline.

Value of α Minimum value of Value of u

h(+ve) for which closed
curve occurs in −ve axis
0.0005 0.0343254079 -0.211397
0.005 0.086987336 -0.30454214
0.05 0.30773621 -0.47066493
0.5 1.7026556091 -0.729614
1.0 3.1119334 -0.80938

Table 1: Minimum values of energy for bound orbits

and the u-nullclines
dVef f
Now from the equation du
= 0 we get:
 
1 −1
u−α 1+ e u =0
u
and solving this equation analytically as well as graphically we get the table below. The
graphical solution is given for a comparison.

11
 u nullcline
1
x -0.0005*exp(-1/x)*(1+1/x)
x -0.005*exp(-1/x)*(1+1/x)
x -0.05*exp(-1/x)*(1+1/x)
x -0.5*exp(-1/x)*(1+1/x)
x -1.0*exp(-1/x)*(1+1/x)

0.5
dV(eff)/du

-0.5

-1
-1 -0.5 0 0.5 1
u

Figure 10: Graphical solution of u-nullclines

Value of α Value of u
0.0005 -0.211397
0.005 -0.30454214
0.05 -0.47066493
0.5 -0.729614
1.0 -0.80938

Table 2:the u-nullclines

vs the coupling constant

This is exactly in agreement with the value provided above in the Table 1.
d2 V
We also note that duef2 f > 0, hence Vef f has minima at those points given in the above
table for different α.
Hence the paths shown in the above figures are stable.

7 Conclusion:
Thus from the above study we have found a set of values of coupling constant for a fixed
energy value above a certain threshold limit. So, it is evident that a proper tuning of the
coupling constant and energy may change the motion of a particle from closed to aperiodic
or vice-versa.
A possible application for this phenomena is the removal of outermost electrons of a heavy
atom. The outermost electrons experience a screened coulomb or yukawa type of interaction
due to the presence of other inner electrons. Thus tuning the strength of external magnetic

12
field, one can tune the angular momentum (J)of the outermost electron which in turn
affects the coupling constant (α = Jmα2 a ). For a fixed energy, a proper choice of magnetic

field can cause easy removal of outermost electrons from their shells. Very heavy elements
or atoms in very excited states, can be exposed to external electric field that will alter the
energy of the outermost electron, keeping the angular momentum conserved, so that its
orbit becomes aperiodic. One can also think of electrons just below the conduction band
in a metal undergoing some periodic orbits, can be taken up to the conduction band by
modulating the coupling constant as mentioned above.
In case of strongly-coupled complex plasma the caging effect on the dust particles can give
rise to some of these closed orbits. This in turn may lead to some oscillating collective
behaviors in complex plasmas.
Another important application of yukawa potential is in astronomy to explain anomaly in
the period of orbits of planets[8], mean motion of planets in solar system[13] and in long
time run of satellites[9]. For astrophysical measurements based on radar signals near the
sun, get also affected by the yukawa correction term in gravitational potential[16]. For all
the above cases the periodicity and closure issues discussed above are quite relevant.

Acknowledgements:
RM and SS are thankful to Mrityunjay Kundu1 , Sayantani Bhattacharyya2 (presently
at IIT Kanpur) and Abhijit Sen1 for their valuable suggestions and discussions. The au-
thors also thank an anonymous referee of Indian Journal of Physics for several insightfull
suggestion.

References:
1. L D Landau and E M Lifshitz Course of Theoretical Physics (Mechanics) (Pergamon
Press : Oxford) Vol 1 Ch 3, Sec 14, p 32 (1969)
2. I Rodriguez and J L Brun Eur. J. Phys. 19 41 (1998).
3. J L Brun and A F Pacheco Celestial Mech. Dyn. Astr 96 311 (2006)
4. M Malek Nonlinear Systems of Ordinary Differential Equations (California State
University, East Bay)
5. S L Ross Differential Equations (John Willey & Sons) Ch 13, Sec 13.3, p 661 (1984)
6. C Sparrow The Qualitative Theory of Ordinary Differential Equations (Lecture Notes
: MA371) p 15 (2009)
7. S H Strogatz Nonlinear Dynamics and Chaos (Levant Books) Ch 8, Sec 8.2, p 253
and Ch 5, Sec 5.2, p 137
8. I Haranas, O Ragos and V Mioc Astrophys Space Sci 332 107 (2011)

13
 9. I Haranas and O Ragos Astrophys Space Sci 331 115 (2011)

10. E Fischbach and C L Talmadge The Search for Non-Newtonian Gravity (Springer) p
113 (1999)

11. J W Moffat and V T Toth ArXiv:0712.1796v5 [gr-qc] (2009)

12. L Iorio Scholarly Research Exchange 2008 Article ID 238385

13. I Haranas, I Kotsireas, G Gomez, M J Fullana, I Gkigkitzis Astrophys Space Sci

361:365 (2016)

14. I Haranas, I Gkigkitzis Astrophys Space Sci 337:693-702 (2012)

15. L Iorio and M L Ruggiero Int. J. Mod. Phys. A 22 (No. 29) 5379-5389 (2007)

16. I Haranas, O Ragos Astrophys Space Sci 334:71-74 (2011)

17. O Bertolami Int. J. Mod. Phys. A 16 (No. 12A) 2003-2012 (2007)

14