# Superradiance and stability of rotating charged black holes in T-duality

Sohan Kumar Jha<EMAIL_ADDRESS>Department of Physics, Chandernagore
College, Chandernagore, Hooghly, West Bengal, India  Kimet Jusufi
<EMAIL_ADDRESS>Physics Department, State University of Tetovo,
Ilinden Street nn, 1200, Tetovo, North Macedonia

###### Abstract

We investigate the shadow images, the relation between Quasinormal Modes
(QNMs) and the shadow radius, and the superradiance effect observed in the
context of a rotating charged black hole under T-duality. Our investigation
places particular emphasis on two key parameters: the electric charge denoted
as $Q$ and the quantum deformed parameter represented by the zero-point
length, $l_{0}$. Our findings reveal a distinct pattern: as the electric
charge increases, the shadow radius experiences a consistent decrease.
Intriguingly, when considering the quantum deformed parameter, we find a
noteworthy phenomenon—a reflecting point. Specifically, we illustrate that the
shadow radius initially increases with an increase in $l_{0}$ and subsequently
decreases. Further analysis involves the computation of eikonal equatorial and
polar QNMs, where a similar reflecting point emerges upon varying $l_{0}$.
This establishes the inverse correlation between QNMs and shadow radius within
our research framework. Our investigation into the effects of Q and $l_{0}$ on
superradiance reveals that the amplification factor initially increases with
$Q$ and $l_{0}$ and then starts decreasing. Moreover, the rotating black holes
in T-duality allows superradiance scattering for a wider range of frequency
than Kerr black holes, making black holes in T-duality brighter than the Kerr
black holes. We also delve into the stability of the combined system of the
rotating black hole and scalar field with the help of black hole bomb
mechanism. It provides a window to observe the impact of parameters Q and
$l_{0}$ on the stability. It shows that the combined system is stable for a
wider regime for the Kerr black hole.

##  I Introduction

Numerous contemporary astrophysical observatories provide strong evidence
supporting the presence of supermassive black holes (SMBHs) at the cores of
giant elliptical and spiral galaxies. One illustrative instance is the
observation indicating the existence of a supermassive black hole with a mass
equivalent to four million solar masses at the heart of our Milky Way galaxy
[1, 2, 3, 4, 5]. These supermassive black holes are characterized by their
masses and spin parameters.

According to Einstein’s theory of relativity, black holes typically contain a
spacetime singularity at their center and an event horizon where gravity is so
intense that nothing, not even electromagnetic radiation or light, can escape.
Due to the external spacetime geometry, a black hole (BH) can capture light
from nearby stars or accretion disks into bound orbits. In simpler terms, each
black hole is characterized by a photon sphere—a collection of light rays
orbiting the BH. It’s noteworthy that the orbit of light is considered
unstable if the photon can either fall into the black hole or escape to
infinity.

Another compelling evidence for the existence of supermassive black hole comes
from the first shadow images at the center of the M87 galaxy that was captured
by the Event Horizon Telescope (EHT) collaboration [6, 7]. Additionally, the
detection of gravitational waves by LIGO contributes to our understanding [8].
Intriguingly, by monitoring the motion of stars in the Galactic Center and
collecting precise measurements of their orbital motion, we can delve into the
nature of black holes and the surrounding spacetime.

Superradiance is a radiation enhancement process where it is possible to
extract energy from the black hole under certain conditions. Penrose first
conceptualized the exciting idea of energy extraction from black holes in
seminal works [9, 10]. Later, Misner, Zel’dovich, and Bekenstein independently
derived the essential condition, $\omega<m\Omega$, for superradiance where
$\omega$, m, and $\Omega$ are the frequency of the incident wave, the
azimuthal quantum number with respect to the rotational axis, and the angular
velocity of the rotating body, respectively [11-16]. The same inequality holds
if we consider the superradiance scattering of electromagnetic and
gravitational waves [17]. It was various investigations in the field of
superradiance that played an instrumental role in the discovery of black hole
evaporation [18]. Studies such as [19, 46] may give an impression of
superradiance’s sole dependence on the horizon boundary condition. However,
numerous studies reveal that the requisite dissipation essential for the
superradiance scattering is provided by the ergoregion [21-24]. It is thus
possible to observe superradiance in the case of horizonless stars [21-27].
Instabilities induced by superradiance in the background spacetime may lead to
a hairy black hole solution, thus providing an excellent window to investigate
the no-hair theorem [28]. Superradiance provides an avenue to probe the
stability of black holes and explore the gravitational collapse in confining
geometries [29].

In this manuscript, our focus is on investigating both the superradiance
effect and the Quasinormal Modes (QNMs) of a charged black hole identified in
T-duality [42]. This particular black hole spacetime is characterized by
parameters such as mass, spin, electric charge, and the quantum-deformed
parameter, also known as the zero-point length parameter. The geometry of this
spacetime is asymptotically flat, the solution is regular at the origin, and
serves as a generalization of the Kerr–Newman solution. In the realm of
astrophysical black holes, the primary influencers are anticipated to be the
mass, spin, and electric charge. However, for smaller black holes and when
exploring the interior structure, quantum effects may play a significant role.
Our exploration delves into understanding these phenomena in the context of
the T-duality framework.

The paper is organized as follows. In Section II, we review the charged black
hole solution in T-duality. In Section III we investigate the shadow images
and the relation between QNMs and shadow radius. In Section IV, we discuss the
superradiance effect and the stability of the black hole. Finally, in Section
V, we comment on our findings.

##  II The BH solution

In this section, we shall review the exact solution for a charged black hole
in T-duality recently found in [42]. In this setup, one starts by assuming the
general solution to be static and spherically symmetric in the form

$ds^{2}=-f(r)\,dt^{2}+\frac{dr^{2}}{g(r)}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})$ (1)

with

$\begin{split}g(r)&=1-\frac{2m(r)}{r}\end{split}$ (2)

where the mass profile is given by

$m(r)=-4\pi\int_{0}^{r}dr^{\prime}(r^{\prime})^{2}\
\mathfrak{T}_{0}^{0}(r^{\prime}).$ (3)

In addition, it is assumed that the effective Einstein field equation holds
[42]

$\displaystyle
R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\kappa\left(\mathfrak{T}_{\mu\nu}^{\mathrm{bare}}+\mathfrak{T}_{\mu\nu}^{\mathrm{em}}\right)$
(4)

where $\kappa=8\pi$, along the modified energy-momentum tensor. Namely we have
a quantum-modified source as well as the modified electromagnetic field

$\mathfrak{T}_{\mu\nu}^{\mathrm{em}}=\frac{1}{4\pi}\left(F_{\mu\sigma}{F_{\nu}}^{\sigma}-\frac{1}{4}g_{\mu\nu}F_{\rho\sigma}F^{\rho\sigma}\right).$
(5)

In T-duality, it has been shown that the gravitational potential is modified
and given by the non-singular expression [42]

$V(r)=-\frac{m_{0}}{\sqrt{r^{2}+l_{0}^{2}}}.$ (6)

From the Newton-Poisson equation, one finds

$\rho_{\mathrm{bare}}(r)=\frac{1}{4\pi}\Delta
V(r)=\frac{3l_{0}^{2}m_{0}}{4\pi{\left(r^{2}+l_{0}^{2}\right)}^{5/2}}.$ (7)

This expression can be viewed as a generalization of the conventional Newton-
Poisson equation, where the source is usually given in terms of the Dirac
delta function. In order to obtain the full solution, one has to consider the
effect of the electromagnetic field using the modified potential [41]

$V_{em}(r)=-\frac{Q}{\sqrt{r^{2}+l_{0}^{2}}},$ (8)

From where one has

$\mathfrak{T}_{0}^{0\
\mathrm{em}}=\frac{Q^{2}r^{2}}{8\pi(r^{2}+l_{0}^{2})^{3}}.$ (9)

As can be seen this expression vanishes at the origin, and it is decaying as
$\sim 1/r^{4}$ at the asymptotic region at infinity. With these results in
hand, one can work out the following exact solution [42]

$\displaystyle f(r)$ $\displaystyle=$ $\displaystyle
g(r)=1-\frac{2m_{0}r^{2}}{\left(r^{2}+l_{0}^{2}\right)^{3/2}}+\frac{5Q^{2}r^{2}}{8(r^{2}+l_{0}^{2})^{2}}$
(10) $\displaystyle+$
$\displaystyle\frac{3Q^{2}l_{0}^{2}}{8(r^{2}+l_{0}^{2})^{2}}-\frac{3Q^{2}}{8l_{0}\,r}\arctan\left(\frac{r}{l_{0}}\right).$

This is a regular spacetime with the ADM mass $M$ given by (see for details
[42])

$M=m_{0}+\frac{3\pi Q^{2}}{32\ l_{0}}.$ (11)

Surprisingly, one finds an intriguing revelation: the mass encompasses an
additional term directly proportional to the regularized self-energy of the
electrostatic field. In the realm of general relativity, this particular term
is typically disregarded, as it is considered integral to the curvature
singularity.

It follows that one can employ the expression (11) to formulate the metric
coefficient in its ultimate representation:

$f(r)=1-\frac{2Mr^{2}}{\left(r^{2}+l_{0}^{2}\right)^{3/2}}+\frac{Q^{2}\
r^{2}}{\left(r^{2}+l_{0}^{2}\right)^{2}}F(r)$ (12)

where

$\displaystyle F(r)$ $\displaystyle=$
$\displaystyle\frac{5}{8}+\frac{3l_{0}^{2}}{8r^{2}}+\frac{3\pi}{16l_{0}}\left(r^{2}+l_{0}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\times$
$\displaystyle\left\\{1-\frac{\left(r^{2}+l_{0}^{2}\right)^{\frac{3}{2}}}{r^{3}}\left[\frac{2}{\pi}\arctan\left(\frac{r}{l_{0}}\right)\right]\right\\}.$

Figure 1: (Left panel) Shadow of the rotating charged black hole for a fixed
$l_{0}=0.1$ (red curve $Q=0.1$, blue curve $Q=0.2$ and black curve $Q=0.3$)
and $M=1$, $a=0.5$ and viewing angle $\theta_{0}=\pi/2$. (Right panel) Shadow
of the rotating charged black hole for a fixed $Q=0.1$ (red curve $l_{0}=0.1$,
blue curve $l_{0}=0.3$ and black curve $Q=0.5$) and $M=1$, $a=0.5$ and viewing
angle $\theta_{0}=\pi/2$.

To get the rotating solution, we can utilize the non-complexification
procedure and we can drop the complexification of the radial coordinate [58].
Finally, in terms of the Kerr-like coordinates, one can show that [42]

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-\frac{\Delta}{\Sigma}(dt-a\sin^{2}\theta
d\phi)^{2}+\frac{\Sigma}{\Delta}dr^{2}+\Sigma\,d\theta^{2}$ (14)
$\displaystyle+$
$\displaystyle\frac{\sin^{2}\theta}{\Sigma}[adt-(r^{2}+a^{2})d\phi]^{2},$

with

$\displaystyle\Delta(r)$ $\displaystyle=$ $\displaystyle
r^{2}-\frac{2Mr^{4}}{\left(r^{2}+l_{0}^{2}\right)^{3/2}}+\frac{Q^{2}r^{4}}{(r^{2}+l_{0}^{2})^{2}}F(r)+a^{2},$

along with

$\displaystyle\Sigma(r,\theta)$ $\displaystyle=$ $\displaystyle
r^{2}+a^{2}\cos^{2}\theta,$ (15)

where we note that $a=L/M$ is the specific angular momentum, and $M$ is the
ADM mass of the black hole.

##  III Black hole shadow and the correspondence with QNMs

In this section, our focus shifts to the study of the evolution of photons in
the vicinity of the spinning black hole in T-duality. We can start the
analysis by using the Hamilton-Jacobi equation

$\displaystyle\frac{\partial
S}{\partial\lambda}=-\frac{1}{2}g^{\mu\nu}\frac{\partial S}{\partial
x^{\mu}}\frac{\partial S}{\partial x^{\nu}},$ (16)

where $\lambda$ is an affine parameter. Next, we can write the Jacobi action
$S$ in the standard form

$S=\frac{1}{2}m_{0}^{2}\lambda-Et+L_{z}\phi+S_{r}(r)+S_{\theta}(\theta),$ (17)

where we shall use the fact that for photons, the rest mass is zero
($m_{0}=0$). Furthermore $E$ and $L_{z}$ are the energy and angular momentum
of the photon, respectively. The functions $S_{r}(r)$ and $S_{\theta}(\theta)$
are introduced, with $S_{r}$ being dependent solely on $r$, and $S_{\theta}$
dependent solely on $\theta$.

Figure 2: (Left panel) The plot of the typical shadow radius with the viewing
angle $\theta_{0}=\pi/2$ as a function of $Q$. We have fixed $l_{0}=0.1$ and
$M=1$, $a=0.5$ and viewing angle $\theta_{0}=\pi/2$. (Right panel) The plot of
the typical shadow radius with the viewing angle $\theta_{0}=\pi/2$ as a
function of $l_{0}$. We have fixed $Q=0.1$, $M=1$, $a=0.5$ and viewing angle
$\theta_{0}=\pi/2$, respectively.

To this end, by substituting the Jacobi action into the Hamilton-Jacobi
equation, we obtain

$\displaystyle S_{r}(r)$ $\displaystyle=$
$\displaystyle\int^{r}\frac{\sqrt{R(r)}}{\Delta}dr,$ (18) $\displaystyle
S_{\theta}(\theta)$ $\displaystyle=$
$\displaystyle\int^{\theta}\sqrt{\Theta(\theta)}d\theta,$ (19)

where

$R(r)=[(r^{2}+a^{2})E-aL_{z}]^{2}-\Delta[\mathcal{K}+(L_{z}-aE)^{2}],$ (20)
$\Theta(\theta)=\mathcal{K}+(a^{2}E^{2}-L^{2}_{z}\csc^{2}\theta)\cos^{2}\theta,$
(21)

where $\mathcal{K}$ is the Carter constant. The Hamilton-Jacobi equation leads
to the derivation of the following four equations of motion governing the
trajectory of photons:

$\displaystyle\rho^{2}\frac{dt}{d\lambda}$ $\displaystyle=$ $\displaystyle
a(L_{z}-aE\sin^{2}\theta)+\frac{r^{2}+a^{2}}{\Delta}[(r^{2}+a^{2})E-aL_{z}],$
$\displaystyle\rho^{2}\frac{d\phi}{d\lambda}$ $\displaystyle=$
$\displaystyle\frac{L_{z}}{\sin^{2}\theta}-aE+\frac{a}{\Delta}[(R^{2}+a^{2})E-aL_{z}],$
$\displaystyle\rho^{2}\frac{dr}{d\lambda}$ $\displaystyle=$
$\displaystyle\sqrt{R(r)},$ $\displaystyle\rho^{2}\frac{d\theta}{d\lambda}$
$\displaystyle=$ $\displaystyle\sqrt{\Theta(\theta)}.$ (22)

Let us define the two impact parameters

$\displaystyle\xi=\frac{L_{z}}{E},\qquad\eta=\frac{\mathcal{K}}{E^{2}},$ (23)

which will be used to investigate the geometric configuration of the black
hole shadow. To achieve this, we employ the unstable condition associated with
circular geodesics given by

$\displaystyle R(r)=0,\qquad\frac{dR(r)}{dr}=0.$ (24)

In the general scenario, photons emitted by a light source undergo deflection
as they traverse near the black hole due to gravitational lensing effects.
Consequently, some photons reach a distant observer after being deflected by
the black hole, while others directly plunge into the black hole. The photons
unable to escape from the black hole collectively shape the observable shadow
of the black hole in the observer’s sky. The boundary of this shadow outlines
the apparent silhouette of the black hole. To examine the characteristics of
the shadow, we adopt celestial coordinates as defined below:

$\displaystyle X$ $\displaystyle=$ $\displaystyle-\xi\csc{\theta_{0}},$ (25)
$\displaystyle Y$ $\displaystyle=$
$\displaystyle\pm\sqrt{\eta+a^{2}\cos^{2}\theta_{0}-\xi^{2}\cot^{2}\theta_{0}}.$
(26)

where [45]

$\displaystyle\xi(r)$ $\displaystyle=$
$\displaystyle\frac{\mathcal{X}_{\text{0}}\Delta^{\prime}_{\text{0}}-2\Delta_{\text{0}}\mathcal{X}^{\prime}_{\text{0}}}{a\Delta^{\prime}_{\text{0}}},$
(27) $\displaystyle\eta(r)$ $\displaystyle=$
$\displaystyle\frac{4a^{2}\mathcal{X}^{\prime
2}_{\text{0}}\Delta_{\text{0}}-\left[\left(\mathcal{X}_{\text{0}}-a^{2}\right)\Delta^{\prime}_{\text{0}}-2\mathcal{X}^{\prime}_{\text{0}}\Delta_{\text{0}}\right]^{2}}{a^{2}\Delta^{\prime
2}_{\text{0}}}.$ (28)

with $\mathcal{X}=r^{2}+a^{2}$. Here, we note that the subscript “0” indicates
the above equations should be evaluated at $r=r_{\text{0}}$, which is a
solution to (24). In Fig. 1, we give the shadow images of the black holes by
varying the electric charge (left panel) and by varying the $l_{0}$ (right
panel). We observe that by keeping a fixed value for $l_{0}$, the shadow
radius decreases with the increase of the electric charge. We can see more
clearly this fact from the plot of the typical shadow radius in Fig. 2 (left
panel). On the other hand, the situation is more interesting when we fix the
electric charge $Q$, and we increase the value of $l_{0}$. As can be seen from
Fig. 2 (right panel), initially, the typical shadow radius increases very
slowly with the increase of $l_{0}$, and then around some critical value
around $l_{0}\sim 0.2$, the shadow radius reaches the maximal value. With a
further increase of $l_{0}$, the shadow radius decreases faster.

###  III.1 The correspondence between shadow radius and eikonal QNMs

In this section, we elaborate the relation between the QNMs, and the shadow
radius (see, [46, 47, 48, 49, 50, 51, 52, 53, 54, 55]). It is well-known that
QNMs modes are associated with the ringdown phase of a black hole, can be
expressed in terms of their real and imaginary parts as
$\omega=\omega_{\Re}-i\omega_{\Im}$. An intriguing observation was shown in
Ref. [51], where it was argued that the QNM frequency of a Kerr black hole in
the eikonal limit is given by:

$\displaystyle\omega_{QNM}=(l+\frac{1}{2})\Omega_{R}-i\gamma_{L}\left(n+\frac{1}{2}\right)$
(29)

with

$\displaystyle\Omega_{R}=\Omega_{\theta}+\frac{m}{l+\frac{1}{2}}\Omega_{prec},$
(30)

where $\Omega_{\theta}$ represents the orbital frequency in the polar
direction, $\Omega_{\text{prec}}$ corresponds to the Lense-Thirring precession
frequency of the orbit plane, $\gamma_{L}$ stands for the Lyapunov exponent of
the orbit, and $n$ denotes the overtone number. In the upcoming discussion, we
would like to study two specific modes, equatorial and polar modes, along with
the effect of electric charge and zero point length.

###  III.2 Viewing angle: $\theta_{0}=\pi/2$

In the general case of a rotating black hole, there is no closed form
expression for the shadow radius. In what follows, we shall consider the case
of a viewing angle: $\theta_{0}=\pi/2$. For simplicity, we shall also consider
equatorial orbit, which can be used to compute the typical shadow radius. We
are also going to use the following fact that the Lense-Thiring precession
frequency is related to the orbital frequency and Keplerian frequency, namely,
the lense–thirring precession frequency for prograde orbits in the limit of a
small perturbation with respect to the equatorial plane is as follows

$\displaystyle\Omega_{prec}=\pm\Omega_{\phi}\mp\Omega_{\theta}$ (31)

where

$\displaystyle\Omega_{\phi}=\frac{-\partial_{r}g_{t\phi}\pm\sqrt{\left(\partial_{r}g_{t\phi}\right)^{2}-(\partial_{r}g_{tt})(\partial_{r}g_{\phi\phi})}}{\partial_{r}g_{\phi\phi}}.$
(32)

In the case of spinning metric in Ref. [48] it was studied the case $m\pm l$,
while in Ref. [52] it was argued that

$\displaystyle\mathcal{K}+L_{z}^{2}\simeq
L^{2}-\frac{a^{2}E^{2}}{2}\left(1-\frac{L_{z}^{2}}{L^{2}}\right).$ (33)

We can rewrite the last equation as follows

$\displaystyle\eta+\xi^{2}\simeq\frac{L^{2}}{E^{2}}-\frac{a^{2}}{2}\left(1-\frac{L_{z}^{2}}{L^{2}}\right).$
(34)

Figure 3: (Left panel) The plot of the real part of the equatorial mode
$\omega_{R}^{+}$ as a function of $Q$. We have fixed $l_{0}=0.1$ and $M=1$,
$a=0.5$ and $l=2$. (Right panel) The plot of the real part of the equatorial
mode $\omega_{R}^{-}$ as a function of $Q$. We have fixed $l_{0}=0.1$, $M=1$,
$a=0.5$ and $l=2$, respectively.

Figure 4: (Left panel) The plot of the real part of the equatorial mode
$\omega_{R}^{+}$ as a function of $l_{0}$. We have fixed $Q=0.1$ and $M=1$,
$a=0.5$ and $l=2$. (Right panel) The plot of the real part of the equatorial
mode $\omega_{R}^{-}$ as a function of $l_{0}$. We have fixed $Q=0.1$, $M=1$,
$a=0.5$ and $l=2$, respectively.

Now, we use the geometric-optics correspondence between the parameters of a
QNMs and the conserved quantities along geodesics. In particular we write:
[51]

$\displaystyle L_{z}$ $\displaystyle\Longleftrightarrow$ $\displaystyle m$
(35) $\displaystyle E$ $\displaystyle\Longleftrightarrow$
$\displaystyle\omega_{\Re}$ (36) $\displaystyle L$
$\displaystyle\Longleftrightarrow$ $\displaystyle l+\frac{1}{2}$ (37)

In the eikonal limit with $m=l\gg 1$, the real part of QNMs is expressed as
$\omega_{\Re}=L\Omega_{R}$. Under this condition, if we introduce
$\mu=m/(l+1/2)=1$, the relation can be established using Eq. (31) with a
positive sign before the orbital frequency:

$\displaystyle\Omega_{prec}=\Omega_{\phi}-\Omega_{\theta},$ (38)

hence one can obtain

$\displaystyle\Omega_{R}=\Omega_{\theta}+\Omega_{prec}=\Omega_{\phi},$ (39)

where the term $\Omega_{\theta}$ cancels out. It means that the real part of
QNMs is related to the Kepler frequency given by

$\displaystyle\omega^{\pm}_{\Re}=(l+\frac{1}{2})\frac{-\partial_{r}g_{t\phi}\pm\sqrt{\left(\partial_{r}g_{t\phi}\right)^{2}-\partial_{r}g_{tt}\partial_{r}g_{\phi\phi}}}{\partial_{r}g_{\phi\phi}}$
(40)

which is valid in the limit $m=l>>1$. Alternatively, we can consider the mode
$m=-l$, which follows that $\mu=-1$, where we can use Eq. (31) with a negative
sign before the Kelperian frequency, and we arrive again at Eq. (40). In this
case, we can use the following definition to specify the shadow radius [56,
43, 44, 48]

$\displaystyle
R_{sh}:=\frac{1}{2}\left(\xi^{+}(r_{0}^{+})-\xi^{-}(r_{0}^{-})\right),$ (41)

provided $\eta(r_{0}^{\pm})=0$. From Eq. (34) one can obtain

$\displaystyle\xi^{\pm}=\pm\sqrt{\frac{(l+\frac{1}{2})^{2}}{\omega^{2}_{\Re}(r_{0}^{\pm})}-\frac{a^{2}}{2}(1-\mu^{2})}.$
(42)

The correspondence is precise if we consider the eikonal limit, that is, if we
set $\mu=\pm 1$ (namely $[(m=\pm l)]$, yielding [43, 44]

$R_{sh}(\mu=\pm
1)=\frac{l+\frac{1}{2}}{2}\left(\frac{1}{\omega_{\Re}(r_{0}^{+})}+\frac{1}{\omega_{\Re}(r_{0}^{-})}\right).$
(43)

The final equation aligns with the expression derived in [48]. Utilizing Eq.
(40) and incorporating the metric functions in the equatorial plane, we can
derive the following expression for the typical shadow radius: [43]

$R_{sh}=\frac{\sqrt{2}}{2}\left(\sqrt{\frac{r_{0}^{+}}{f^{\prime}(r)|_{r_{0}^{+}}}}+\sqrt{\frac{r_{0}^{-}}{f^{\prime}(r)|_{r_{0}^{-}}}}\right).$
(44)

The preceding equation corresponds precisely to the result established in Ref.
[48], where the points $r_{0}^{\pm}$ were determined by solving the equation:

$r_{0}^{2}-\frac{2r_{0}}{f^{\prime}(r)|_{r_{0}^{\pm}}}f(r_{0})\mp
2a\sqrt{\frac{2r_{0}}{f^{\prime}(r)|_{r_{0}^{\pm}}}}=0.$ (45)

As we elaborated above, in the eikonal limit, we can use Eq. (40), which can
be further simplified as follows

$\displaystyle\omega_{\Re}^{\pm}=(l+\frac{1}{2})\frac{1}{a\pm\sqrt{\frac{2r_{0}^{\pm}}{f^{\prime}(r)|_{r_{0}^{\pm}}}}}.$
(46)

Using the last equation, in Fig. 3, we plotted our result for the real part of
equatorial QNMs. As we can see, with the increase of charge, the value of
$\omega_{\Re}$ increases. This confirms the inverse relation between QNMs and
shadow radius, where we found that shadow radius decreases with the increase
of the electric charge (see Fig. 2). Furthermore, when varying $l_{0}$, we
observe a reflecting point as initially, $\omega_{\Re}$ decreases with the
increase of $l_{0}$, reaching some minimal value and then increasing faster by
increasing $l_{0}$ (see Fig. 4, right panel). To conclude, it was shown that
the reflecting point was resent in the case of shadow radius and QNMs; hence,
the correspondence between QNMs and shadow radius holds.

###  III.3 Viewing angle: $\theta_{0}=0$ & $\theta_{0}=\pi$

We will now focus on the polar orbit, in that case we have to use $\theta=0$.
In this specific scenario, it is natural to calculate the shadow radius with
viewing angles for the observer, denoted as $\theta_{0}=0$ and
$\theta_{0}=\pi$. Utilizing celestial coordinates, we can derive:

$\displaystyle X^{2}+Y^{2}=\eta+\xi^{2}+a^{2}\cos^{2}\theta_{0}.$ (47)

This shows that, when considering a viewing angle $\theta_{0}=0$ (or
$\theta_{0}=\pi$), the shadow retains circular disk shape. Therefore, it is
natural to employ the following definition [56]

$\displaystyle
R_{sh}:=\sqrt{a^{2}+\eta(r_{0})},\,\,\,\text{with}\,\,\xi(r_{0})=0.$ (48)

In the context of the polar orbit, the azimuthal angular momentum is zero,
i.e., $L_{z}=0$. Under these conditions, and considering the requirements for
the existence of circular geodesics, specifically $\dot{r}^{2}$, we deduce
from the radial geodesics:

$(r^{2}+a^{2})^{2}-[r^{2}f(r)+a^{2}]R_{sh}^{2}=0,$ (49)

and

$4r(r^{2}+a^{2})-2rf(r)R_{sh}^{2}-r^{2}f^{\prime}(r)R_{sh}^{2}=0,$ (50)

where the definition $R_{sh}^{2}=\mathcal{K}/E^{2}+a^{2}$ (see in particular
[55]) has been used. By utilizing Eq. (49), the typical shadow radius can be
defined as $R_{sh}\equiv(R_{sh}^{+}-R_{sh}^{-})/2$. Consequently, we obtain:

$R_{sh}=\frac{a^{2}+r^{2}}{\sqrt{r^{2}f(r)+a^{2}}}|_{r_{0}},$ (51)

In this context, $r_{0}$ can be determined by solving the algebraic equation:

$\displaystyle(a^{2}+r_{0}^{2})^{2}-\frac{4[r_{0}^{2}f(r_{0})+a^{2}](a^{2}+r_{0}^{2})}{r_{0}f^{\prime}(r_{0})+2f(r_{0})}=0.$
(52)

On the other hand, using Eq. (34) and the geometric-optic correspondence, we
get the relation between the shadow radius and the real part of QNMs in terms
of

$\displaystyle
R_{sh}=\sqrt{\frac{(l+1/2)^{2}}{\omega^{2}_{\Re}(r_{0})}+\frac{a^{2}}{2}}.$
(53)

Figure 5: (Left panel) The plot of the real part of the polar mode
$\omega_{R}$ as a function of $Q$. We have fixed $l_{0}=0.1$ and $M=1$,
$a=0.5$ and $l=2$. (Right panel) The plot of the real part of the polar mode
$\omega_{R}$ as a function of $l_{0}$. We have fixed $Q=0.1$ and $M=1$, along
with $a=0.5$ and $l=2$.

Finally, using Eq. (51) along with the last equation, we get for the real part
of QNMs a compact equation

$\omega_{\Re}=\left(l+\frac{1}{2}\right)\sqrt{\frac{2(r^{2}f(r)+a^{2})\zeta}{\zeta}},$
(54)

where $\zeta=2r^{4}+4a^{2}r^{2}+a^{4}-a^{2}r^{2}f(r)$. In Fig. 5 (left panel),
we plotted our result for the real part of polar QNMs. As we can see, with the
increase of charge, the value of $\omega_{R}$ increases. On the other hand,
for the polar QNMs, when varying $l_{0}$ we observe a reflecting point, that
is, initially $\omega_{\Re}$ decreases with the increase of $l_{0}$, reaching
some minimal value and then increasing faster by increasing $l_{0}$ (see Fig.
5, right panel). It should be noted that for larger values of $l$, we should
obtain more accurate values for QNMs, as the correspondence is accurate in the
eikonal limit.

##  IV Superradiance

In this section, we intend to study the superradiance scattering of a massive
scalar field $\Phi$ off a rotating black hole in T-duality. To this end, the
Klein-Gordon equation in curved spacetime is taken into consideration. It is
given by

$\displaystyle\left(\bigtriangledown_{\alpha}\bigtriangledown^{\alpha}+\mu^{2}\right)\Phi(t,r,\theta,\phi)=0,$

which can be also written as

$\displaystyle\left[\frac{-1}{\sqrt{-g}}\partial_{\sigma}\left(g^{\sigma\tau}\sqrt{-g}\partial_{\tau}\right)+\mu^{2}\right]\Phi(t,r,\theta,\phi)=0.$

Here, $\mu$ represents the mass of the scalar field $\Phi$. With the help of
the standard separation of variables method and using the following ansatz
with the standard Boyer-Lindquist coordinates $(t,r,\theta,\phi)$,

$\displaystyle\Phi(t,r,\theta,\phi)=R_{\omega jm}(r)\Theta(\theta)e^{-i\omega
t}e^{im\phi},$ (55)

with $\quad j\geq 0,\quad-j\leq m\leq j,\quad\omega>0$. We will be able to
separate Eq. (IV) into radial and angular parts. In Eq. (55), $R_{\omega
jm}(r)$ and $\Theta(\theta)$ are the radial function and oblate spheroidal
wave function, respectively. The symbols $j$, $m$, and $\omega$ represent the
angular eigenfunction, angular quantum number, and the positive frequency of
the field under investigation as measured by a far-away observer,
respectively. The Eq. (IV), with the help of the ansatz (55), yields following
radial and angular equations:

$\displaystyle\frac{d}{dr}(\Delta\frac{dR_{\omega
jm}(r)}{dr})+(\frac{((r^{2}+a^{2}\omega-am)^{2}}{\Delta})R_{\omega lm}(r)$
$\displaystyle-(\mu^{2}r^{2}+j(j+1)+a^{2}\omega^{2}-2m\omega a)R_{\omega
lm}(r)=0,$ (56)
$\displaystyle\sin\theta\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta_{\omega
jm}(\theta)}{d\theta}\right)+$ (57) $\displaystyle+$
$\displaystyle\left(j(j+1)\sin^{2}\theta-\left(\left(a\omega\sin^{2}\theta-m\right)^{2}\right)\right)\Theta_{\omega
jm}(\theta)$ $\displaystyle+$ $\displaystyle
a^{2}\mu^{2}\sin^{2}\theta\cos^{2}\theta~{}\Theta_{\omega jm}(\theta)=0.$

###  IV.1 Amplification factor for superradiance

Since our intention is to investigate the scattering of $\Phi$, we will only
consider Eq. (56) for the rest of the paper. Following articles [30, 31], the
general solution of the radial equation (56) can be obtained. Next, We
introduce a Regge-Wheeler-like coordinate $r_{*}$ defined by

$\displaystyle r_{*}\equiv\int dr\frac{r^{2}+a^{2}}{\Delta},$ (58)

with $r_{*}\rightarrow-\infty$ at event horizon and $r_{*}\rightarrow\infty$
at infinity. We then define a new radial function $\mathcal{P}_{\omega
jm}\left(r_{*}\right)=\sqrt{r^{2}+a^{2}}R_{\omega jm}(r)$ to transform the
radial equation (56) into the desired shape. A few steps of algebra provide
the following equation:

$\frac{d^{2}\mathcal{P}_{\omega lm}\left(r_{*}\right)}{dr_{*}^{2}}+V_{\omega
jm}(r)\mathcal{P}_{\omega jm}\left(r_{*}\right)=0,$ (59)

where $V_{\omega jm}(r)$ is the effective potential given by

$\displaystyle V_{\omega jm}(r)$ $\displaystyle=$
$\displaystyle\left(\omega-\frac{ma}{r^{2}+a^{2}}\right)^{2}-\frac{\Delta}{\left(r^{2}+a^{2}\right)^{2}}\Big{[}j(j+1)$
(60) $\displaystyle+$ $\displaystyle a^{2}\omega^{2}-2ma\omega+\mu^{2}r^{2}$
$\displaystyle+$
$\displaystyle\sqrt{r^{2}+a^{2}}\frac{d}{dr}\left(\frac{r\Delta}{\left(r^{2}+a^{2}\right)^{\frac{3}{2}}}\right)\Big{]},$

So, the resulting effect is that we now study the scattering of the scalar
field $\Phi$ under the effective potential (60). The asymptotic behavior of
the potential (60) at the event horizon and at spatial infinity play a crucial
role in obtaining the asymptotic solutions of the radial equation (56). The
limiting value of the potential at the event horizon is

$\displaystyle\lim_{r\rightarrow r_{eh}}V_{\omega
jm}(r)=\left(\omega-m\Omega_{t}\right)^{2}\equiv
k_{eh}^{2},\quad\Omega_{t}=\frac{a}{r_{eh}^{2}+a^{2}},$ (61)

and the limiting value at spatial infinity is

$\lim_{r\rightarrow\infty}V_{\omega jm}(r)=\omega^{2}-\mu^{2}\equiv
k_{\infty}^{2}.$ (62)

The potential has distinct, constant limiting values at the extreme points.
With the asymptotic values of the potential at our disposal, we are now able
to obtain the asymptotic radial solutions. After a few steps of algebra, we
obtain

$R_{\omega
jm}(r)\rightarrow\left\\{\begin{array}[]{cl}\frac{\mathcal{N}_{in}^{eh}e^{-ik_{eh}r_{*}}}{\sqrt{r_{eh}^{2}+a^{2}}}&\text{
for }r\rightarrow r_{eh}\\\
\mathcal{N}_{in}^{\infty}\frac{e^{-ik_{\infty}r_{*}}}{r}+\mathcal{N}_{ref}^{\infty}\frac{e^{ik_{\infty}r_{*}}}{r}&\text{
for }r\rightarrow\infty\end{array}\right.$ (63)

Here, $\mathcal{N}_{in}^{eh}$ and $\mathcal{N}_{in}^{\infty}$ correspond to
the amplitude of the incoming scalar wave at the event horizon($r_{eh}$) and
infinity $(\infty)$, respectively. $\mathcal{N}_{ref}^{\infty}$ represents the
amplitude of the reflected part of the scalar wave at infinity $(\infty)$.
Matching Wronskinas at the event horizon and infinity, we obtain

$\left|\mathcal{N}_{ref}^{\infty}\right|^{2}=\left|\mathcal{N}_{in}^{\infty}\right|^{2}-\frac{k_{eh}}{k_{\infty}}\left|\mathcal{N}_{in}^{eh}\right|^{2}.$
(64)

Figure 6: Variation of percentage amplification factors against $\omega$ for
different values of $l_{0}$ with $Q=0.4$.

The above relation shows that the superradiance scattering of the scalar field
is possible when $frac{k_{eh}}{k_{\infty}}<0$, i.e., $\omega<m\Omega_{T}$. To
obtain the analytical expressions of amplitudes, we employ the asymptotic
matching method proposed in [32, 33]. In this procedure, we obtain separate
solutions for the two overlapping regions: one is the near region where
$r-r_{eh}\ll\omega$ and another is the far region where $r-r_{eh}>>1$. We then
expand the resulting solutions to match them in the overlapping region which
finally yields the desired expressions of amplitudes. Following the procedure
and assumptions given in [34, 35, 36], we get the amplification factor as

$Z_{jm}\equiv\frac{\left|\mathcal{N}_{ref}^{\infty}\right|^{2}}{\left|\mathcal{N}_{in}^{\infty}\right|^{2}}-1,$
(65)

where

$\displaystyle\mathcal{N}_{in}^{\infty}$ $\displaystyle=$
$\displaystyle\frac{b(-2i)^{-\frac{1+\sqrt{1+4j(j+1)}}{2}}}{\sqrt{(\omega^{2}-\mu^{2})}}\cdot\frac{\Gamma(\sqrt{1+4j(j+1)})\Gamma(1+\sqrt{1+4j(j+1)})}{\Gamma\left(\frac{1+\sqrt{1+4j(j+1)}}{2}-2iB\right)\left(\Gamma\left(\frac{1+\sqrt{1+4j(j+1)}}{2}\right)\right)^{2}}\times$
$\displaystyle\Gamma(1-2iB)k_{t}^{\frac{1-\sqrt{1+4j(j+1)}}{2}}+\frac{b(-2i)^{\frac{\sqrt{1+4j(j+1)}-1}{2}}}{\sqrt{(\omega^{2}-\mu^{2})}}\times$
$\displaystyle\frac{\Gamma(1-\sqrt{1+4j(j+1)})\Gamma(-\sqrt{1+4j(j+1)})}{\left(\Gamma\left(\frac{1-\sqrt{1+4j(j+1)}}{2}\right)\right)^{2}\Gamma\left(\frac{1-\sqrt{1+4j(j+1)}}{2}-2iB\right)}\Gamma(1-2iB)k_{t}^{\frac{1+\sqrt{1+4j(j+1)}}{2}},$

and

$\displaystyle\mathcal{N}_{ref}^{\infty}$ $\displaystyle=$
$\displaystyle\frac{b(2i)^{-\frac{1+\sqrt{1+4j(j+1)}}{2}}}{\sqrt{(\omega^{2}-\mu^{2})}}\cdot\frac{\Gamma(\sqrt{1+4j(j+1)})\Gamma(1+\sqrt{1+4j(j+1)})}{\Gamma\left(\frac{1+\sqrt{1+4j(j+1)}}{2}-2iB\right)\left(\Gamma\left(\frac{1+\sqrt{1+4j(j+1)}}{2}\right)\right)^{2}}\times$
$\displaystyle\Gamma(1-2iB)k_{t}^{\frac{1-\sqrt{1+4j(j+1)}}{2}}+\frac{b(2i)^{\frac{\sqrt{1+4j(j+1)}-1}{2}}}{\sqrt{(\omega^{2}-\mu^{2})}}\times$
$\displaystyle\frac{\Gamma(1-\sqrt{1+4j(j+1)})\Gamma(-\sqrt{1+4j(j+1)})}{\left(\Gamma\left(\frac{1-\sqrt{1+4j(j+1)}}{2}\right)\right)^{2}\Gamma\left(\frac{1-\sqrt{1+4j(j+1)}}{2}-2iB\right)}\Gamma(1-2iB)k_{t}^{\frac{1+\sqrt{1+4j(j+1)}}{2}}.$

Here, $B=\frac{\omega-m\Omega_{t}}{r_{eh}-r_{ch}}r_{eh}^{2}$ and
$k_{t}=\left(r_{eh}-r_{ch}\right)\sqrt{\omega^{2}-\mu^{2}}$, $r_{ch}$ being
the Cauchy horizon. As evident from the Eq. (65), we have superradiance when
$Z_{jm}>0$. Since cases $m<0$ are non-superradiant, we ignore them here as we
are interested in the occurrence of superradiance. We intend to study the
effect of the ”Zero-point length” $l_{0}$ and the charge Q. We have taken
$M=1$, $\mu=0.1$, and $a=0.5$.

Figure 7: Variation of percentage amplification factors against $\omega$ for
different values of $Q$ with $l_{0}=0.1$.

Figure 8: Comparison between Kerr black hole and the rotating black hole in
T-duality.

Through the display of plots (6, 7), we illustrate the effect of parameters
$l_{0}$ and Q on the superradiance scattering. Both the plots show that the
amplification factor initially increases with an increase in $l_{0}$ and Q and
then starts decreasing. This trend is also visible in Fig. (8), where we make
a comparison between the Kerr black hole and the rotating black hole in
T-duality. It is also observed that the black hole in T-duality allows wider
range of frequency to be amplified. For those frequencies that lie between the
critical frequencies of the Kerr black hole and black hole in T-duality, waves
get absorbed by the Kerr black hole but superradiantly scattered by the black
hole in T-duality. This has larger impact on the luminosity. To illustrate
this fact further, we consider a monochromatic massive scalar field. The
outgoing energy flux of the field, measured by an observer at infinity, is
given by [37]

$\dot{E}=\frac{\omega
k_{\infty}}{2}\left|\mathcal{N}_{ref}^{\infty}\right|^{2}=\frac{\omega
k_{\infty}}{2}\frac{\left|\mathcal{N}_{ref}^{\infty}\right|^{2}}{\left|\mathcal{N}_{in}^{\infty}\right|^{2}}\left|\mathcal{N}_{in}^{\infty}\right|^{2}$
(68)

We take the ratio of energy fluxes for the Kerr black hole and the black hole
in T-duality at the respective critical superradiant frequency to make a
meaningful comparison. The ratio of critical energy fluxes for massless scalar
fields with same $\left|\mathcal{N}_{in}^{\infty}\right|$ is

$\frac{\dot{E}^{crit}_{T}}{\dot{E}^{crit}_{Kerr}}=\frac{\Omega_{t}^{2}}{\Omega_{kerr}^{2}}>1$
(69)

The ratio being greater than unity shows that the black hole in T-duality
would be brighter than the Kerr black hole. To illustrate the dependency of
the ratio on the parameters $l_{0}$ and Q, we plot the ratio as a function of
Q for different values of $l_{0}$ in Fig. (9). Here, the spin parameter is
fixed at $0.5$.

Figure 9: Ratio of the outgoing energy fluxes as a function of charge Q for
different values of $l_{0}$.

###  IV.2 Superradiant instability for the rotating black hole in T-duality

In this section, we investigate the superradiant instability of the combined
system of rotating black hole and scalar field due to the phenomenon known as
black hole bomb [38, 39]. In this phenomenon, the superradiant field gets
trapped by a potential well and gets repeatedly amplified, resulting in
exponential growth of the amplitude and hence superradiance instability.

From equation (56), we have

$\displaystyle\Delta\frac{d}{dr}\left(\Delta\frac{dR_{\omega
jm}}{dr}\right)+\mathcal{K}R_{\omega jm}=0,$ (70)

where for a slowly rotating black hole $(\hat{a}\omega\ll 1)$

$\mathcal{K}\equiv\left(\left(r^{2}+a^{2}\right)\omega-
ma\right)^{2}+\Delta\left(2ma\omega-j(j+1)-\mu^{2}r^{2}\right).$

Following the black hole bomb mechanism, solutions for the radial equation
(70) is given by

$R_{\omega
jm}\sim\left\\{\begin{array}[]{ll}e^{-i(\omega-m\Omega_{t})r_{*}}&\text{ as
}r\rightarrow r_{eh}\left(r_{*}\rightarrow-\infty\right)\\\
\frac{e^{-\sqrt{\mu^{2}-\omega^{2}}r_{*}}}{r}&\text{ as
}r\rightarrow\infty\left(r_{*}\rightarrow\infty\right)\end{array}\right.$

We now introduce a new radial function:

$\psi_{\omega jm}\equiv\sqrt{\Delta}R_{\omega jm},$

that converts the radial equation (70) into

$\left(\frac{d^{2}}{dr^{2}}+\omega^{2}-V\right)\psi_{\omega jm}=0.$

with

Figure 10: Parameter space($m\Omega$-$\mu$) for massive scalar field where
colored area represents region with stable dynamics and non-colored area
represents region with unstable dynamics. Here $Q=0.4$.

Figure 11: Parameter space($m\Omega$-$\mu$) for massive scalar field where
colored area represents region with stable dynamics and non-colored area
represents region with unstable dynamics. Here $l_{0}=0.1$.

$\omega^{2}-V=\frac{\mathcal{K}+\frac{1}{4}(\frac{d\Delta}{dr})^{2}-\frac{\Delta}{2}\frac{d^{2}\Delta}{dr^{2}}}{\Delta^{2}},$

which is the Regge-Wheel equation. Retaining terms up to the order
$\mathcal{O}\left(1/r\right)$, we obtain the asymptotic form of the effective
potential $V(r)$ as

$V(r)=\mu^{2}-\frac{4M\omega^{2}}{r}+\frac{2M\mu^{2}}{r}.$

We impose the condition $V^{\prime}\rightarrow 0^{+}$ as $r\rightarrow\infty$
essential for the realization of black hole bomb mechanism [40]. Taking into
consideration the fact that the superradiance amplification of scattered waves
occurs when $\omega<m\Omega_{t}$, we obtain

$\frac{\mu}{\sqrt{2}}<\omega^{2}<m\Omega_{t},$

in which the integrated system of the rotating black hole in T-duality and the
massive scalar field may experience a superradiant instability, known as the
black hole bomb. The dynamics of the massive scalar field in black hole will
remain stable when $\mu\geq\sqrt{2}m\Omega_{t}$. Plots (10, 11) show that the
parameter space $(m\Omega-\mu)$ over which the combined system of black hole
and scalar field remains stable decreases with an increase in either $l_{0}$
or Q.

##  V Final remarks

In this current study, we have explored the shadow images, the correlation
between QNMs and the shadow radius, as well as the investigation of the
superradiance effect and the stability of charged black holes within the
context of T-duality. Our findings indicate that the shadow radius decreases
as the electric charge increases. Particularly intriguing is the observation
regarding the quantum deformed parameter $l_{0}$, where a distinct phenomenon
emerges—a reflecting point is identified. At this point, the shadow radius
initially increases with an increase in $l_{0}$ and subsequently decreases.
Furthermore, our exploration extended to the analysis of equatorial and polar
QNMs, revealing a notable pattern: the real part of QNMs amplifies with an
increase in charge. Nevertheless, an intriguing turning point is uncovered
when varying $l_{0}$. Consequently, this shows the inverse relationship
between QNMs and shadow radius in our investigation. Broadly speaking, in the
context of astrophysical black holes, the predominant influence is attributed
to the electric charge, while the role of $l_{0}$ is expected to be important
only for smaller black holes. It it still an open question if such quantum
effects can have an impact for large black holes.

For the superradiance the amplification factor shows the existence of
reflecting points for both the parameters. The critical frequency, however,
increases with an increase in the electric charge and the quantum deformed
parameter. This has a direct bearing on the luminosity of the black hole as
observed at infinity. A wider range of allowed frequency for superradiance
makes black holes in T-duality brighter than the Kerr black hole.

## References

  * [1] The EHT Collaboration, Astrophys. J. 875, L1 (2019).
  * [2] The EHT Collaboration, Astrophys. J. Lett. 910, L13 (2021).
  * [3] K. Akiyama et al. [Event Horizon Telescope], Astrophys.J.Lett. 930 (2022) L12
  * [4] K. Akiyama et al. [Event Horizon Telescope], Astrophys.J.Lett. 930 (2022) L13
  * [5] K. Akiyama et al. [Event Horizon Telescope], Astrophys.J.Lett. 930 (2022) L14
  * [6] K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875, L1 (2019)
  * [7] K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875, no.1, L2 (2019)
  * [8] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016)
  * [9] R. Penrose, Gen. Relativ. Gravit. 34, 1141 (2002).
  * [10] R. Penrose and R. M. Floyd, Nature (London) 229, 177 (1971).
  * [11] C.W. Misner, Phys. Rev. Lett. 22, 1071 (1969).
  * [12] Ya. B. Zel’dovich, JETP Lett. 14, 180 (1971).
  * [13] Ya. B. Zel’dovich, J. Exp. Teor. Phys. 35, 1085 (1972),
  * [14] S.W. Hawking, Phys. Rev. Lett. 26, 1344 (1971).
  * [15] J. D. Bekenstein, Phys. Rev. D 7, 949 (1973).
  * [16] J. D. Bekenstein and M. Schiffer, Phys. Rev. D 58, 064014 (1998).
  * [17] S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * [18] S.W. Hawking, Commun. Math. Phys. 43, 199 (1975); 46,206(E) (1976).
  * [19] M. Richartz, S. Weinfurtner, A. J. Penner, and W. G. Unruh, Phys. Rev. D 80, 124016 (2009).
  * [20] V. Cardoso and P. Pani, Classical Quantum Gravity 30,045011 (2013).
  * [21] G. Eskin, Rev. Math. Phys. 28, 1650025 (2016).
  * [22] R. Vicente, V. Cardoso, and J. C. Lopes, Phys. Rev. D 97,084032 (2018).
  * [23] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, 1(2015).
  * [24] D.N. Page, Phys. Rev. D 13 (1976) 198.
  * [25] M. Richartz and A. Saa, Phys. Rev. D 88, 044008(2013).
  * [26] V. Cardoso, R. Brito, and J. L. Rosa, Phys. Rev. D 91,124026 (2015).
  * [27] K. Glampedakis, S. J. Kapadia, and D. Kennefick, Phys.Rev. D 89, 024007 (2014).
  * [28] A. Rahmani, M. Khodadi, M. Honardoost, and H. R. Sepangi, Nucl. Phys. B960, 115185 (2020).
  * [29] Cardoso V., Dias O. J., Phys.Rev. D70 (2004) 084011,
  * [30] V. B. Bezerra , H. S. Vieira and A. A. Costa: Class.Quantum Grav. 31 045003 (2014)
  * [31] G. V. Kraniotis: Class.Quant.Grav. 33 225011 (2016)
  * [32] A. A. Starobinsky, Sov. Phys. JETP 37, 28 (1973), http://www.jetp.ras.ru/cgi-bin/e/index/e/37/1/p28?a=list.
  * [33] A. A. Starobinsky and S.M. Churilov, Zh. Eksp. Teor. Fiz. 65, 3 (1973), http://www.jetp.ras.ru/cgibin/ e/index/e/38/1/p1?a=list [Sov. Phys. JETP 38, 1 (1973)].
  * [34] S. Kr. Jha, A. Rahaman Eur. Phys. J. C (2022) 82:728
  * [35] S. Kr. Jha et. al. PHYSICAL REVIEW D 107, 084052 (2023)
  * [36] S. Kr. Jha, A. Rahaman Physics of the Dark Universe 42 (2023) 101327
  * [37] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, 1 (2015).
  * [38] S.A. Teukolsky, W.H. Press, Astrophys. J. 193 (1974) 443.
  * [39] W.H. Press, S.A. Teukolsky, Nature 238 (1972) 211–212.
  * [40] S. Hod, Phys. Lett. B 708 (2012) 320 [arXiv:1205.1872 [gr-qc]].
  * [41] P. Gaete and P. Nicolini, Phys. Lett. B 829 (2022), 137100
  * [42] P. Gaete, K. Jusufi and P. Nicolini, Phys. Lett. B 835 (2022), 137546
  * [43] K. Jusufi, Phys. Rev. D 101 (2020) no.12, 124063
  * [44] K. Jusufi, M. Azreg-Aïnou, M. Jamil and Q. Wu, Universe 8 (2022) no.4, 210
  * [45] R. Shaikh, Phys. Rev. D 100, 024028 (2019).
  * [46] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79, 064016 (2009).
  * [47] I. Z. Stefanov, S. S. Yazadjiev and G. G. Gyulchev, Phys. Rev. Lett. 104, 251103 (2010).
  * [48] K. Jusufi, Phys. Rev. D 101, 124063 (2020).
  * [49] K. Jusufi, Phys. Rev. D 101 (2020) no.8, 084055
  * [50] R. A. Konoplya and Z. Stuchlík, Phys. Lett. B 771 (2017), 597-602
  * [51] H. Yang, D. A. Nichols, F. Zhang, A. Zimmerman, Z. Zhang and Y. Chen, Phys. Rev. D 86, 104006 (2012).
  * [52] H. Yang, Phys. Rev. D 103, 084010 (2021).
  * [53] B. Cuadros-Melgar, R. D. B. Fontana and J. de Oliveira, Phys. Lett. B 811, 135966 (2020).
  * [54] B. Mashhoon, Phys. Rev. D 31, 290 (1985).
  * [55] S. R. Dolan, Phys. Rev. D 82, 104003 (2020).
  * [56] X. H. Feng and H. Lu, arXiv:1911.12368 [gr-qc].
  * [57] M. Zhang and M. Guo, Eur. Phys. J. C 80, 790 (2020).
  * [58] M. Azreg-Aïnou, Phys. Rev. D 90 (2014) no.6, 064041


11institutetext: Centre for Nonlinear Science and Engineering, School of
Electrical and Electronics Engineering, SASTRA Deemed University,
Thanjavur-613401, Tamil Nadu, India 22institutetext: Centre for Nonlinear
Systems, Chennai Institute of Technology, Chennai-600069, Tamil Nadu, India
33institutetext: Centre for Nonlinear Science and Engineering, School of
Electrical and Electronics Engineering, SASTRA University, Thanjavur-613401,
Tamil Nadu, India 44institutetext: Department of Nonlinear Dynamics, School of
Physics, Bharathidasan University, Tiruchirappalli - 620024, Tamil Nadu, India
55institutetext: Department of Nonlinear Dynamics, School of Physics,
Bharathidasan University, Tiruchirappalli - 620024, Tamil Nadu, India

# A Class of Isochronous and Non-Isochronous Nonlinear Oscillators

J. Ramya Parkavi _E-mail:_<EMAIL_ADDRESS>11 R. Mohanasubha _E-mail:_
<EMAIL_ADDRESS>V. K. Chandrasekar _E-mail:_
<EMAIL_ADDRESS>(Corresponding Author)33 M. Senthilvelan _E-mail:_
<EMAIL_ADDRESS>M. Lakshmanan _E-mail:_<EMAIL_ADDRESS>
(Received: date / Revised version: date)

###### Abstract

In this work, we present a method of generating a class of nonlinear ordinary
differential equations (ODEs), representing the dynamics of appropriate
nonlinear oscillators, that have the characteristics of either amplitude
independent frequency of oscillations or amplitude dependent frequency of
oscillations from the integrals of the simple harmonic oscillator equation. To
achieve this, we consider the case where the integrals are in the same form
both for the linear and the nonlinear oscillators in either of the cases. We
also discuss the method of deriving the associated integrals and the general
solution in harmonic form for both the types. We demonstrate the applicability
of this method up to $2N$ coupled first order nonlinear ODEs in both the
cases. Further, we illustrate the theory with an example in each case.

##  1 Introduction

In recent years, a great deal of interest has been shown in identifying
nonlinear dynamical systems that posses isochronous property RIac -bud . The
systems whose frequencies of oscillations are independent of the amplitudes
similar to that of a linear harmonic oscillator are being termed as
isochronous systems. Isochronous systems also occur in several physical
situations, for example, starting from the linear harmonic oscillator to
nonlinear equations like Li$\acute{e}$nard type equations. Isochronicity
phenomenon has been widely studied not only for its impact in stability
theory, but also for its relationship with bifurcation and boundary value
problems. However, only a class of nonlinear systems exhibits isochronous
oscillations. To identify and study the behaviour of such isochronous
nonlinear dynamical systems, several methods have been developed in the
literature, see for example Refs. cho1 -guha . In this work, we present a
novel method of generating nonlinear ordinary differential equations (ODEs)
that have the characteristics of isochronous behaviour. We derive explicit
general solutions and present integrals of the identified isochronous
nonlinear systems. This class of isochronous nonlinear oscillators can be
related to linear harmonic oscillator by introducing a nonlocal transformation
in the dependent variable.

Several nonlinear dynamical systems exhibit non-isochronous oscillations. In
these oscillators one can see that the frequency of oscillations depends on
the amplitude. For example, the anharmonic oscillator equation
$\ddot{x}+c_{1}x+c_{2}x^{3}=0$, with $c_{1}>0$, $c_{2}>0$, admits the solution
$x(t)=Acn(\omega t,k)$, $\omega=\sqrt{c_{1}+c_{2}A^{2}}$, where $cn$ is the
Jacobi elliptic function of modulus $k=\sqrt{c_{2}A^{2}/2(c_{1}+c_{2}A^{2})}$.
While deriving this solution we have considered the initial conditions
$x(0)=A$, $\dot{x}(0)=0$ mlsr . Unlike the linear harmonic oscillator equation
the angular frequency $\omega$ is now not a constant but depends on the
amplitude $A$. These non-isochronous dynamical systems exhibit a variety of
dynamics differing from the isochronous systems. However, a class of non-
isochronous dynamical systems also admit periodic oscillations of harmonic
type. A notable example in this category is the Mathews-Lakshmanan nonlinear
oscillator equation, $(1+\lambda x^{2})\ddot{x}-\lambda
x\dot{x}^{2}+\omega_{0}^{2}x=0$, where $\lambda$ and $\omega_{0}$ are
constants laksh ; laksh75 . This nonlinear/nonpolynomial oscillator equation
has the solution $x(t)=A\cos\Omega t$, $\Omega=\omega_{0}/\sqrt{1+\lambda
A^{2}}$. We have obtained this solution for the initial conditions $x(0)=A$,
$\dot{x}(0)=0$ mlsr . The angular frequency ($\Omega$) is amplitude dependent
and hence the Mathews-Lakshmanan oscillator is a non-isochronous system.
Because of the amplitude dependent frequency property, various responses have
also been exhibited for different input conditions by the underlying dynamical
system. This kind of nonisochornicity nature can be seen in various physical
and biological systems. In particular, one finds that strong nonisochronicity
produces various synchronization patterns as well as long irregular transient
phase dynamics in networks of such oscillators in the case where the long
transient dynamics has connection with the functioning of biological systems
last1 ; last2 ; last3 .

Upon observing that certain non-isochronous systems can also have periodic
oscillations of harmonic type as in the case of the above Mathews-Lakshmanan
oscillator (and different from that of the anharmonic oscillator mentioned
earlier), we pose the following question. Is it possible to generate nonlinear
isochronous and non-isochronous systems whose solutions are periodic functions
of harmonic type through a single method? Our following investigation shows
that one can indeed generate a class of isochronous and non-isochronous
nonlinear oscillators which admits periodic solutions of harmonic type from
the known integrals of the simple harmonic oscillator equation.

By investigating the time series plot of the harmonic oscillator and nonlinear
isochronous oscillators we find that both of them exhibit the same frequency
of oscillations but they differ only in their amplitudes. Since the angular
frequency originates from the second integral of the harmonic oscillator one
may consider the case where the second integral of the nonlinear systems may
also be in the same form as that of the linear harmonic oscillator equation.
Starting from this observation we have identified a class of nonlinear ODEs
that exhibits isochronous oscillations from the second integral of the linear
harmonic oscillator equation.

On the other hand, in non-isochronous systems there is an interdependency
between the amplitude and the frequency of oscillations. Interestingly in the
case of non-isochronous systems one can observe that the first integral has
the same structure as that of the harmonic oscillator. The only difference we
come across in the integrals of nonlinear systems is that the presence of an
additional function in terms of the phase. This function has the variables $x$
and $\dot{x}$. By carefully examining this function we find that it is nothing
but the Hamiltonian. Interestingly while deriving the general solution from
the obtained integrals the same function appears as amplitude with a square
root. This in turn establishes an interrelation between the amplitude and the
frequency of oscillations. In this way we generate non-isochronous systems
from the known integrals of the linear harmonic oscillator equation.

As far as our knowledge goes, the generalized version of the Mathews-
Lakshmanan oscillator equation to $2N$ dimensions is considered for the first
time. Another notable point is that the dynamics of the $2N$-coupled Mathews-
Lakshmanan oscillator is studied for the non-identical system parameters (i.e.
$\omega_{1}\neq\omega_{2}$). Further, we have presented the general solution
of the $2N$-coupled Mathews-Lakshmanan oscillator equation. Also we have
studied the integrability of the modified Emden equation with more generalized
parameters which has also not been studied previously in the literature.

We organize our work as follows. In Sec. 2, we present the method of
identifying nonlinear oscillator equations which exhibit isochronous and non-
isochronous properties. In Sec. 3, we consider the isochronous case and
discuss the method of deriving the solution for a system of two coupled first
order ODEs. We also extend the method of finding the solutions and the
integrals to $2N$ coupled first-order ODEs and demonstrate it with suitable
examples. In Sec. 4, we consider the non-isochronous case and present the
method of finding the solution for a system of two coupled first order ODEs.
We prove the method with a physically important example, namely Mathews-
Lakshmanan oscillator. Then we extend the method to $2N$-coupled first order
ODEs and present the general solution for the $2N$-coupled Mathews-Lakshmanan
oscillator equation. Finally, we present our conclusions in Sec. 5.

##  2 Generating nonlinear systems having harmonic and anharmonic behaviour

Let us consider a system of two coupled first order nonlinear ODE’s of the
form

$\displaystyle\dot{x}=P(x,y),~{}~{}\dot{y}=Q(x,y),$ (1)

where $P$ and $Q$ are functions of $x$ and $y$ and the overdot denotes
differentiation with respect to $t$. We assume that the above system of
nonlinear ODEs can be connected to the linear harmonic oscillator equation,
that is

$\displaystyle\dot{u}=\omega v,~{}\dot{v}=-\omega u,$ (2)

where $\omega$ is the angular frequency of oscillation. Equation (2) admits an
integral of the form

$I_{1l}=-\omega t+\arctan\Big{(}\frac{v}{u}\Big{)}.$ (3)

Suppose the nonlinear ODE (1) also admits the same form of integral, that is

$\displaystyle I_{1nl}$ $\displaystyle=$ $\displaystyle-\omega
t+\arctan(\frac{y}{x}).$ (4)

Equations (3) and (4) imply that ${v}/{u}={y}/{x}$. To satisfy this equality,
we correlate the variables in the form

$x=\frac{u}{h(x,y)},~{}y=\frac{v}{h(x,y)},$ (5)

where $h(x,y)$ is an arbitrary function of $x$ and $y$.

Differentiating Eq. (5) with respect to $t$ and substituting back (2) and (5)
in the resultant equations and simplifying them, we obtain

$\displaystyle P(x,y)=\omega y-f(x,y)x,~{}~{}~{}Q(x,y)=-\omega x-f(x,y)y,$ (6)

where $f(x,y)=\frac{\dot{h}}{h}$.

The linear harmonic oscillator Eq.(2) has another integral (energy integral or
first integral) of the form $I_{2l}=u^{2}+v^{2}$. Let the nonlinear and linear
ODEs $(\ref{eq:1})$ and $(\ref{eq:2})$ have the same form of the first
integral. By equating these two integrals, we find

$\displaystyle u^{2}+v^{2}=x^{2}+y^{2}.$ (7)

Differentiating (7) with respect to $t$ and replacing $\dot{x}=P(x,y)$ and
$\dot{y}=Q(x,y)$, we obtain

$xP(x,y)+yQ(x,y)=0.$ (8)

To fulfill the condition (8), we choose the following specific forms for $P$
and $Q$, namely

$\displaystyle P(x,y)=\omega f(x,y)y,~{}~{}~{}Q(x,y)=-\omega f(x,y)x,$ (9)

where $f$ is an arbitrary function of $x$ and $y$. We note that the form
considered for $P$ and $Q$ in (9) is only a specific choice. One may consider
other forms for $P$ and $Q$ also. However, in this work, we confine our
attention only to this choice.

The nonlinear ODE (1) which has the right hand side of the form (6) admits
isochronous behaviour whereas the nonlinear ODEs (1) which have the right hand
side of the form (9) exhibit non-isochronous behaviour, as we see below in
sections $3$ and $4$.

##  3 Isochronous systems

###  3.1 2-coupled isochronous systems

To begin with, we focus our attention on the class of systems described by the
form (6), that is

$\displaystyle\dot{x}=\omega y-{f(x,y)}x,~{}\dot{y}=-\omega x-{f(x,y)}y.$ (10)

From (10), we can deduce the following expression, that is

$\displaystyle y\dot{x}-x\dot{y}=\omega{(x^{2}+y^{2})}.$ (11)

Equation (11) can be rewritten in the form

$\displaystyle\frac{1}{\omega}\frac{d}{dt}\bigg{(}\frac{x}{y}\bigg{)}=1+\bigg{(}{\frac{x}{y}\bigg{)}^{2}}.$
(12)

By introducing a new variable $\theta$ as

$\frac{x}{y}=\tan\theta,$ (13)

Eq. (12) can be reduced to

$\displaystyle\dot{\theta}=\omega\Rightarrow\theta=\omega t+\delta,$ (14)

where $\delta$ is an integration constant. From (14) it is clear that for a
class of ODEs considered in (10), the angular velocity $(\dot{\theta})$ or the
angular frequency $(\omega)$ is nothing but a constant.

Even though the angular frequency is a constant, depending on the form of
$f(x,y)$, the amplitude of oscillations may decay or grow and therefore one
has to identify a suitable form of $f(x,y)$ which allows periodic
oscillations. After analyzing the outcome, we find that whenever the function
$f(x,y)$ satisfies the condition

$\displaystyle\int_{0}^{2\pi/\omega}f(\sin(\omega
t^{\prime}+\delta),\cos(\omega t^{\prime}+\delta))dt^{\prime}=0$ (15)

then the underlying system (10) admits an isochronous center. In other words,
whenever the function $f(x,y)$ satisfies the above condition, the system (10)
exhibits periodic oscillations and the associated dynamical system is an
isochronous one. Differing from the above, the forms which do not satisfy the
condition $(\ref{mod10aa})$ display other types of oscillatory behaviour.

###  3.2 $2N$ coupled isochronous systems

The methodology given above can be generalized to higher order ODEs also. To
demonstrate this, let us consider a system of $2N$-coupled first order
nonlinear ODEs of the form [2]

$\displaystyle\dot{x}_{i}$ $\displaystyle=$
$\displaystyle\omega_{i}y_{i}-f(\overline{x},\overline{y})x_{i},$
$\displaystyle\dot{y}_{i}$ $\displaystyle=$
$\displaystyle-\omega_{i}x_{i}-f(\overline{x},\overline{y})y_{i},\hskip
28.45274pti=1,2,...,N,$ (16)

where
$f(\overline{x},\overline{y})=f(x_{1},x_{2},...,x_{N},y_{1},y_{2},...,y_{N})$.
We note that while extending the analysis to higher order ODEs we restrict the
function $f(\bar{x},\bar{y})$ which appears in (16) to be of the same form.
Multiplying the first and second equations of (16) by $y_{i}$ and $x_{i}$
respectively, and subtracting the resulting equations suitably, we obtain

$\displaystyle
y_{i}\dot{x}_{i}-x_{i}\dot{y}_{i}=\omega_{i}(x^{2}_{i}+y^{2}_{i}),~{}~{}~{}i=1,2,...,N.$
(17)

Dividing Eqs. (17) by $y^{2}_{i}$, $i=1,...,N,$ and rewriting it, we find

$\displaystyle\frac{d}{dt}\Big{(}\frac{x_{i}}{\omega_{i}y_{i}}\Big{)}=1+\Big{(}\frac{x_{i}}{y_{i}}\Big{)}^{2},\hskip
28.45274pti=1,2,...,N.$ (18)

Upon introducing the angle variable
$\displaystyle\theta_{i}=\tan^{-1}\Big{(}\frac{x_{i}}{y_{i}}\Big{)}$, Eq. (18)
becomes

$\displaystyle\dot{\theta}_{i}=\omega_{i},\hskip 28.45274pti=1,2,...,N.$ (19)

Equation (19) reveals that the angular frequency is independent of the
amplitudes of oscillations, irrespective of the form of
$f(\overline{x},\overline{y})$.

In the following, we demonstrate the existence of $2N$ independent integrals
for the system of equations (16). To construct the integrals of Eq. (16), we
recall

$\displaystyle\frac{y_{i}}{x_{i}}=\cot(\omega_{i}t+\delta_{i}),\hskip
28.45274pti=1,2,...,N.$ (20)

From Eq. (20), we can identify $N$ time dependent integrals which are of the
form

$\displaystyle\delta_{i}=\cot^{-1}\bigg{[}\frac{y_{i}}{x_{i}}\bigg{]}-\omega_{i}t,\hskip
28.45274pti=1,2,...,N.$ (21)

To determine the remaining $N$ integrals of Eq. (16), we proceed as follows.
From Eq.(16), we can identify

$\displaystyle\frac{x_{j}\dot{x}_{j}+y_{j}\dot{y}_{j}}{x_{j}^{2}+y_{j}^{2}}=\frac{x_{N}\dot{x}_{N}+y\dot{y}_{N}}{x_{N}^{2}+y_{N}^{2}}.\hskip
28.45274ptj=1,2,...,N-1.$ (22)

Rewriting Eq. (22) as
$\frac{d}{dt}\left[\log(x_{j}^{2}+y_{j}^{2})/2\right]=\frac{d}{dt}\left[\log(x_{N}^{2}+y_{N}^{2})/2\right]$
and integrating and rewriting the resultant integrals for $j=1,2,\cdots,N-1$,
we obtain $N-1$ time independent integrals of the form

$\displaystyle I_{j}=\frac{x_{j}^{2}+y_{j}^{2}}{x_{N}^{2}+y_{N}^{2}}.\hskip
28.45274ptj=1,2,...,N-1.$ (23)

By utilizing the obtained $(N-1)$ integrals we can reduce the set of $2N$
first order equations to a system of two coupled first order ODEs in the
variables $x_{N}$ and $y_{N}$. In the following, we demonstrate this.

From Eq.(20), we find

$\frac{x_{i}}{y_{i}}=\tan(\theta_{i}),\;\theta_{i}=\omega_{i}t+\delta_{i}.$
(24)

Using the well-known trigonometric identities, that is
$\sin{\theta_{i}}=\frac{\tan{\theta_{i}}}{\sqrt{1+\tan^{2}{\theta_{i}}}}$ and
$\cos{\theta_{i}}=\frac{\cot{\theta_{i}}}{\sqrt{1+\cot^{2}{\theta_{i}}}}$, we
can identify

$\displaystyle\sin(\omega_{i}t+\delta_{i})=\frac{x_{i}}{\sqrt{x^{2}_{i}+y^{2}_{i}}},$
(25)
$\displaystyle\cos(\omega_{i}t+\delta_{i})=\frac{y_{i}}{\sqrt{x^{2}_{i}+y^{2}_{i}}},\hskip
28.45274pti=1,2,...,N.$ (26)

From (25), we find

$\displaystyle x_{i}(t)$ $\displaystyle=$
$\displaystyle\sin(\omega_{i}t+\delta_{i})(\sqrt{x^{2}_{i}+y^{2}_{i}})$ (27)
$\displaystyle=$
$\displaystyle\sqrt{I_{i}}\sin(\omega_{i}t+\delta_{i})(\sqrt{x^{2}_{N}+y^{2}_{N}})$
$\displaystyle=$
$\displaystyle\sqrt{I_{i}}\frac{\sin(\omega_{i}t+\delta_{i})}{\sin(\omega_{N}t+\delta_{N})}x_{N}.$

Similarly

$\displaystyle y_{i}(t)$ $\displaystyle=$
$\displaystyle\sqrt{I_{i}}\frac{\cos(\omega_{i}t+\delta_{i})}{\cos(\omega_{N}t+\delta_{N})}y_{N}.$
(28)

With the help of Eqs. (27) and (28) one can replace the variables $x_{i}$ and
$y_{i}$, $i=1,2,\dots N-1$, in terms of $t,x_{N}$ and $y_{N}$. As a result, we
find

$\displaystyle f(\overline{x},\overline{y})=f(x_{1},x_{2},\dots
x_{N},y_{1},y_{2},\dots y_{N})=f(t,x_{N},y_{N}).$ (29)

Now confining our attention on the last two ODEs in $(\ref{mod43})$, we find

$\displaystyle\dot{x}_{N}$ $\displaystyle=$
$\displaystyle\omega_{N}y_{N}-f(t,x_{N},y_{N})x_{N},$ (30)
$\displaystyle\dot{y}_{N}$ $\displaystyle=$
$\displaystyle-\omega_{N}x_{N}-f(t,x_{N},y_{N})y_{N}.$ (31)

The above equation cannot be integrated explicitly for an arbitrary form of
$f(t,x_{N},y_{N})$. For particular forms of $f(t,x_{N},y_{N})$, we can solve
the above equation. The theory developed above helps to identify a class of
nonlinear isochronous systems and their general solution.

###  3.3 Examples

In the following, with appropriate choices of the function $f(x,y)$, we
demonstrate through specific examples how specific nonlinear integrable
dynamical systems corresponding to $N=1,2$ and arbitrary integer can be
identified which admit isochronous oscillations. The crucial property which we
impose on the choice of the function $f(x,y)$ is that it should satisfy the
integral condition given by Eq.(15). The condition will also hold good when
$N$ is arbitrary, as for example in Eqs.(30) and (31), the condition (15) will
be imposed on the form of $f(t,x_{N},y_{N})$, as illustrated in the examples
given below.

####  3.3.1 Example 1

To start with we consider the simplest case, namely $N=2$. Based on our
general prescription given above, we consider the choice

$\displaystyle{f(x,y)}=a_{1}x^{q}+a_{2}y^{q},$ (32)

where $q$ is a positive integer and $a_{1}$ and $a_{2}$ are constants.
Substituting the expression $x=y\tan\theta$, with $\theta=\omega t+\delta$,
the second expression in Eq. (10) can be brought to the form

$\displaystyle\dot{y}={(a_{1}\tan^{q}{(\omega
t+\delta)}+a_{2})y^{q+1}}-\omega\tan{(\omega t+\delta)}y.$ (33)

Equation (33) is of Bernoulli form. The general solution of Eq. (33) reads

$\displaystyle y(t)=\frac{\cos{(\omega
t+\delta)}}{{\Big{(}I_{1}-q\int\big{(}a_{1}{\sin}^{q}{(\omega
t+\delta)}+a_{2}{\cos}^{q}{(\omega t+\delta)\big{)}dt\Big{)}^{\frac{1}{q}}}}}$
(34)

from which the expression for $x(t)$ can be deduced as

$\displaystyle x(t)=\frac{\sin{(\omega
t+\delta)}}{{\Big{(}I_{1}-q\int\big{(}a_{1}{\sin}^{q}{(\omega
t+\delta)}+a_{2}{\cos}^{q}{(\omega
t+\delta)\big{)}dt\Big{)}^{\frac{1}{q}}}}}.$ (35)

The integral which is appearing in the denominator in expressions Eqs. (34)
and (35) can be evaluated depending on whether $q$ is an odd or even integer.

We note here that the function given in (32) satisfies the condition (15) only
for odd values of $q$ so that the solution exhibits periodic behaviour. For
even value of $q$ the function (32) does not satisfy the condition (15) and
the solution exhibits damped behaviour. In the following, we present the
solutions for both the cases.

####  3.3.2 Case 1: $q$ is an odd integer $(q=2m+1)$

For this choice, we obtain the solution $x(t)$ and $y(t)$ in the form
tablebook

$\displaystyle x(t)=\sin{(\omega
t+\delta)}{\Bigg{[}I_{1}-q(b_{1}(t)+b_{2}(t))\Bigg{]}}^{\frac{-1}{2m+1}},$
$\displaystyle y(t)=\cos{(\omega
t+\delta)}{\Bigg{[}I_{1}-q(b_{1}(t)+b_{2}(t))\Bigg{]}}^{\frac{-1}{2m+1}},$
(36)

with

$\displaystyle
b_{1}(t)=\frac{a_{1}}{\omega}\bigg{[}\frac{-\cos\tau}{2m+1}\bigg{(}\sin^{2m}\tau+\sum_{k=0}^{m-1}s_{k}^{1}\sin^{2m-2k-2}\tau\bigg{)}\bigg{]},$
$\displaystyle
b_{2}(t)=\frac{a_{2}}{\omega}\bigg{[}\frac{\sin\tau}{2m+1}\bigg{(}\cos^{2m}\tau+\sum_{k=0}^{m-1}s_{k}^{1}\cos^{2m-2k-2}\tau\bigg{)}\bigg{]}.$
(37)

Here $\tau=\omega t+\delta$ and
$s_{k}^{1}=(2^{k+1}m(m-1)\ldots(m-k))/((2m-1)(2m-3)\ldots(2m-2k-1))$.

####  3.3.3 Case 2: $q$ is even integer $(q=2m)$

In this case, the solution (33) takes the form tablebook

$\displaystyle x(t)=\sin{(\omega
t+\delta)}{\Bigg{[}I_{1}-q(b_{3}(t)+b_{4}(t))\Bigg{]}}^{\frac{-1}{2m}},$
$\displaystyle y(t)=\cos{(\omega
t+\delta)}{\Bigg{[}I_{1}-q(b_{3}(t)+b_{4}(t))\Bigg{]}}^{\frac{1}{2m}},$ (38)

where

$\displaystyle
b_{3}(t)=\frac{a_{1}}{\omega}\bigg{[}\frac{-\cos\tau}{2m}\bigg{(}\sin^{2m-1}\tau+\sum_{k=0}^{m-1}s_{k}^{2}\sin^{2m-2k-1}\tau\bigg{)}+\frac{(2m-1)!!}{2^{m}m!}\tau\bigg{]},$
$\displaystyle
b_{4}(t)=\frac{a_{2}}{\omega}\bigg{[}\frac{\sin\tau}{2m}\bigg{(}\cos^{2m-1}\tau+\sum_{k=0}^{m-1}s_{k}^{2}\cos^{2m-2k-1}\tau\bigg{)}+\frac{(2m-1)!!}{2^{m}m!}\tau\bigg{]}.$
(39)

Here $s_{k}^{2}=((2m-1)(2m-3)\ldots(2m-2k+1))/(2^{k}(m-1)(m-2)\ldots(m-k))$.

For odd and even values of $q$ the system (10) with the form of $f(x,y)$ given
in (32) exhibits periodic and damped oscillations, respectively. The solution
plot given in Figure 1 confirms the same.

Figure 1: Time series plot of equation (10) (black line) for (a) $q=1$, (b)
$q=2$, (c) $q=3$, and (d) $q=4$ show both periodic (for $q$ odd values) and
damped (for $q$ even values) oscillation. Red/grey line is the time series
plot of equation (2). Initial conditions and parameter values are
$x(0)=u(0)=-0.5$, $y(0)=v(0)=0.5$ and $a=b=\omega=1.0$.

Case: $f(x,y)=ax+by,~{}q=1$

For the choice $f(x,y)=ax+by$ ($a_{1}=a$, $a_{2}=b$ and $q=1$ in Eq. (32)),
the coupled first order nonlinear ODEs (1) read as

$\displaystyle\dot{x}=\omega y-(ax+by)x,~{}~{}~{}\dot{y}=-\omega x-(ax+by)y.$
(40)

The above system is a special case of the most general quadratic system, that
was studied by Loud loud .

Rewriting Eq. (40) in second order form, we find

$\displaystyle\ddot{x}+3ax\dot{x}+a^{2}x^{3}+\omega x(\omega-
bx)+\frac{1}{(\omega-bx)}(2b\dot{x}^{2}+3abx^{2}\dot{x}+a^{2}bx^{4})=0.$ (41)

To the best of our knowledge, the general solution of the above generalized
version of the modified Emden equation is reported here for the first time.
The general solution of Eq. (41) is given by (vide Eqs.(34) and (35))

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle\frac{\omega\sin(\omega
t+\delta)}{\omega I_{1}+a\cos(\omega t+\delta)-b\sin(\omega t+\delta)},$
$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle\frac{\omega\cos(\omega
t+\delta)}{\omega I_{1}+a\cos(\omega t+\delta)-b\sin(\omega
t+\delta)},~{}~{}\mathrm{I_{1}:constant}.$ (42)

By substituting $f(x,y)=x$ in (6), we find $\dot{x}=\omega y-x^{2}$,
$\dot{y}=-\omega x-xy$. Rewriting this as a single second order ODE in the
variable $x$, we obtain the modified Emden equation annapre , that is
$\ddot{x}+3x\dot{x}+x^{3}+\omega^{2}x=0$. The solution of this equation can be
extracted from the Eqs. $(\ref{18})$ and (35) by fixing $q=1$ and $a_{2}=0$.
The resultant form agrees with the one reported in the literature annapre .

For the choice $q=2$ in Eq. (32), the resultant system turns out to be a cubic
system. The isochronous cases of these systems have already been studied in
the literature, see for example Refs. chav ; chou2 .

Figure 2: Time series plot of equations (2) and (40) (red and black,
respectively) exhibiting periodic oscillations with the same frequency with
different amplitudes for the initial conditions (a) $x(0)=u(0)=-0.5$ and
$y(0)=v(0)=0.5$ and (b) $x(0)=u(0)=-1$ and $y(0)=v(0)=0.5$. Corresponding
phase space portraits are given in (b) and (d) for $a=b=\omega=1.0$.

Figure 2 shows the periodic oscillations admitted by Eqs. (2) and (40) for two
different sets of initial conditions. One may note that the frequency of
oscillations of (41) is independent of the initial conditions as in the case
of the linear harmonic oscillator. However, the amplitude of oscillations of
the nonlinear system differs from the linear one.

####  3.3.4 Example 2

For the case $n=2$, the equation of motion can be written in the form

$\displaystyle\dot{x}_{1}$ $\displaystyle=$
$\displaystyle\omega_{1}y_{1}-f(\bar{x},\bar{y})x_{1},$
$\displaystyle\dot{x}_{2}$ $\displaystyle=$
$\displaystyle\omega_{2}y_{2}-f(\bar{x},\bar{y})x_{2},$
$\displaystyle\dot{y}_{1}$ $\displaystyle=$
$\displaystyle-\omega_{1}x_{1}-f(\bar{x},\bar{y})y_{1},$
$\displaystyle\dot{y}_{2}$ $\displaystyle=$
$\displaystyle-\omega_{2}x_{2}-f(\bar{x},\bar{y})y_{2},$ (43)

where $f(\bar{x},\bar{y})=f(x_{1},x_{2},y_{1},y_{2})$. To proceed further, we
assume a specific form for the function $f(x,y)$, namely

$f=a_{1}x_{1}^{q}+b_{1}y_{1}^{q}+a_{2}x_{2}^{q}+b_{2}y_{2}^{q},$ (44)

where $q$ is an integer and $a_{i}$’s, $b_{i}$’s, $i=1,2$, are constants. Now
the equation of motion for $q=1$ turns out to be

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle
w_{1}y_{1}-(a_{1}x_{1}+b_{1}y_{1}+a_{2}x_{2}+b_{2}y_{2})x_{1},$
$\displaystyle\dot{y}_{1}$ $\displaystyle=$ $\displaystyle-
w_{1}x_{1}-(a_{1}x_{1}+b_{1}y_{1}+a_{2}x_{2}+b_{2}y_{2})y_{1},$
$\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle
w_{2}y_{2}-(a_{1}x_{1}+b_{1}y_{1}+a_{2}x_{2}+b_{2}y_{2})x_{2},$
$\displaystyle\dot{y}_{2}$ $\displaystyle=$ $\displaystyle-
w_{2}x_{2}-(a_{1}x_{1}+b_{1}y_{1}+a_{2}x_{2}+b_{2}y_{2})y_{2},$ (45)

where $w_{i}$ and $a_{j}$’s, $b_{j}$’s $i=1,2$ and $j=1,2,3,4$, are all real
arbitrary constants. Following the procedure discussed in the previous section
with $q=1$, we can identify the general solution of $(\ref{52})$ in the form

$\displaystyle
x_{1}(t)=\frac{\sqrt{I_{1}}w_{1}w_{2}\sin\theta_{1}}{G_{3}},~{}~{}y_{1}(t)=\frac{\sqrt{I_{1}}w_{1}w_{2}\cos\theta_{1}}{G_{3}},$
(46) $\displaystyle
x_{2}(t)=\frac{w_{1}w_{2}\sin\theta_{2}}{G_{3}},~{}~{}y_{2}(t)=\frac{w_{1}w_{2}\cos\theta_{2}}{G_{3}},$
(47)

where $\theta_{1}=w_{1}t+\delta_{1}$, $\theta_{2}=w_{2}t+\delta_{2}$,
$G_{3}=I_{2}w_{1}w_{2}-a_{1}\sqrt{I_{1}}w_{2}\cos\theta_{1}-a_{2}w_{1}\cos\theta_{2}+b_{1}\sqrt{I_{1}}w_{2}\sin\theta_{1}+b_{2}w_{1}\sin\theta_{2}$.

For the choice, $q=2$, we consider an equation which is of the form

$\displaystyle\dot{x}_{1}=w_{1}y_{1}-(a_{1}x_{1}^{2}+b_{1}y_{1}^{2}+a_{2}x_{2}^{2}+b_{2}y_{2}^{2})x_{1},$
$\displaystyle\dot{y}_{1}=-w_{1}x_{1}-(a_{1}x_{1}^{2}+b_{1}y_{1}^{2}+a_{2}x_{2}^{2}+b_{2}y_{2}^{2})y_{1},$
$\displaystyle\dot{x}_{2}=w_{2}y_{2}-(a_{1}x_{1}^{2}+b_{1}y_{1}^{2}+a_{2}x_{2}^{2}+b_{2}y_{2}^{2})x_{2},$
$\displaystyle\dot{y}_{2}=-w_{2}x_{2}-(a_{1}x_{1}^{2}+b_{1}y_{1}^{2}+a_{2}x_{2}^{2}+b_{2}y_{2}^{2})y_{2},$
(48)

where $w_{i}$ and $a_{j}$’s, $b_{j}$’s $i=1,2$ and $j=1,2,3,4$, are all real
arbitrary constants. The general solution of $(\ref{neweq123})$ can be
identified in a straightforward manner as

$\displaystyle
x_{1}(t)=\frac{\sqrt{2I_{1}w_{1}w_{2}}\sin\theta_{1}}{G_{4}},~{}~{}y_{1}(t)=\frac{\sqrt{2I_{1}w_{1}w_{2}}\cos\theta_{1}}{G_{4}},$
$\displaystyle
x_{2}(t)=\frac{\sqrt{2w_{1}w_{2}}\sin\theta_{2}}{G_{4}},~{}~{}y_{2}(t)=\frac{\sqrt{2w_{1}w_{2}}\cos\theta_{2}}{G_{4}},$
(49)

where $\theta_{1}=w_{1}t+\delta_{1}$, $\theta_{2}=w_{2}t+\delta_{2}$,
$G_{4}=\big{[}2w_{1}w_{2}I_{2}+2((a_{1}+b_{1})I_{1}w_{2}\theta_{1}+(a_{2}+b_{2})w_{1}\theta_{2})-(a_{1}-b_{1})I_{1}w_{2}\sin
2\theta_{1}+(b_{2}-a_{2})w_{1}\sin 2\theta_{2}\big{]}^{\frac{1}{2}}$, $I_{1}$
and $I_{2}$ are the integration constants. The above solution leads to damped
type solution which is similar to the Figure 1 for the values $q=2$ and $q=4$.

####  3.3.5 Example 3

For $2N$ coupled system, we consider a specific form for $f$ as

$\displaystyle{f(\overline{x},\overline{y})}=\sum_{j=1}^{N}a_{j}x_{j}^{q}+b_{j}y_{j}^{q},$
(50)

where $a_{j}$’s and $b_{j}$’s are constants. The equation of motion becomes

$\displaystyle\dot{x}_{i}$ $\displaystyle=$
$\displaystyle\omega_{i}y_{i}-\big{[}\sum_{j=1}^{N}a_{j}x_{j}^{q}+b_{j}y_{j}^{q}\big{]}x_{i},$
$\displaystyle\dot{y}_{i}$ $\displaystyle=$
$\displaystyle-\omega_{i}x_{i}-\big{[}\sum_{j=1}^{N}a_{j}x_{j}^{q}+b_{j}y_{j}^{q}\big{]}y_{i},\hskip
28.45274pti=1,2,...,N,$ (51)

For this form of $f$, we have to find $x_{N}$ and $y_{N}$ to write the general
solution of (51) (vide Eqs.(27) and (28)). Eq.(31) becomes

$\displaystyle\dot{y}_{N}(t)$ $\displaystyle=$
$\displaystyle-\omega_{N}\tan(\omega_{N}t+\delta_{N})y_{N}-\sum_{i=1}^{N}\bigg{(}a_{j}I_{j}^{\frac{q}{2}}\frac{\sin^{q}(\omega_{i}t+\delta_{i})}{\cos^{q}(\omega_{N}t+\delta_{N})}$
(52)
$\displaystyle+b_{j}I_{j}^{\frac{q}{2}}\frac{\cos^{q}(\omega_{i}t+\delta_{i})}{\cos^{q}(\omega_{N}t+\delta_{N})}\bigg{)}y_{N}^{q+1}.$

Equation $(\ref{mod43ab})$ is of Bernoulli type and the general solution is
given by

$\displaystyle y_{N}(t)$ $\displaystyle=$
$\displaystyle\frac{\cos(\omega_{N}t+\delta_{N})}{\bigg{[}I_{N}-q\int\sum_{j=1}^{N}(a_{j}I_{j}^{\frac{q}{2}}\sin^{q}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{\frac{q}{2}}\cos^{q}(\omega_{i}t+\delta_{i}))dt\bigg{]}^{\frac{1}{q}}}.$
(53)

Substituting this back in $x_{N}=y_{N}\tan(\omega_{N}+\delta_{N})$, we find
the following expression for $x_{N}(t)$, that is

$\displaystyle x_{N}(t)$ $\displaystyle=$
$\displaystyle\frac{\sin(\omega_{N}t+\delta_{N})}{\bigg{[}I_{N}-q\int\sum_{j=1}^{N}(a_{j}I_{j}^{\frac{q}{2}}\sin^{q}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{\frac{q}{2}}\cos^{q}(\omega_{i}t+\delta_{i}))dt\bigg{]}^{\frac{1}{q}}}.$
(54)

Here also, the integral which is appearing in the denominator can be evaluated
depending on whether $q$ is an odd or even integer. By substituting the above
form of $x_{N}$ and $y_{N}$ in Eqs. (27) and (28), we get the general solution
of Eq. (51). Now the general solution turns out to be

$\displaystyle x_{i}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{I_{i}}\sin(\omega_{i}t+\delta_{i})}{\bigg{[}I_{N}-q\int\sum_{j=1}^{N}(a_{j}I_{j}^{\frac{q}{2}}\sin^{q}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{\frac{q}{2}}\cos^{q}(\omega_{i}t+\delta_{i}))dt\bigg{]}^{\frac{1}{q}}},$
$\displaystyle y_{i}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{I_{i}}\cos(\omega_{i}t+\delta_{i})}{\bigg{[}I_{N}-q\int\sum_{j=1}^{N}(a_{j}I_{j}^{\frac{q}{2}}\sin^{q}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{\frac{q}{2}}\cos^{q}(\omega_{i}t+\delta_{i}))dt\bigg{]}^{\frac{1}{q}}}.$
(55)

For odd integers ($q=2m+1$), the integral appearing in the denominator can be
integrated to yield

$\displaystyle\int\sum_{j=1}^{N}(a_{j}I_{j}^{\frac{2m+1}{2}}\sin^{2m+1}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{\frac{2m+1}{2}}\cos^{2m+1}(\omega_{i}t+\delta_{i}))dt$
$\displaystyle=\sum_{j=1}^{N}\frac{a_{j}}{\omega_{j}}I_{j}^{\frac{2m+1}{2}}\bigg{[}\frac{-\cos\tau_{j}}{2m+1}\bigg{(}\sin^{2m}\tau_{j}+\sum_{k=0}^{m-1}s_{k}^{1}\sin^{2m-2k-2}\tau_{j}\bigg{)}\bigg{]}$
$\displaystyle+\sum_{j=1}^{N}\frac{b_{j}}{\omega_{j}}I_{j}^{\frac{2m+1}{2}}\bigg{[}\frac{\sin\tau_{j}}{2m+1}\bigg{(}\cos^{2m}\tau_{j}+\sum_{k=0}^{m-1}s_{k}^{1}\cos^{2m-2k-2}\tau_{j}\bigg{)}\bigg{]},$
(56)

where $\tau_{j}=\omega_{j}t+\delta_{j}$ and
$s_{k}^{1}=(2^{k+1}m(m-1)\ldots(m-k))/((2m-1)(2m-3)\ldots(2m-2k-1))$.

For even integers ($q=2m$), we find

$\displaystyle\int\sum_{j=1}^{N}(a_{j}I_{j}^{m}\sin^{2m}(\omega_{i}t+\delta_{i})+b_{j}I_{j}^{m}\cos^{2m}(\omega_{i}t+\delta_{i}))dt$
$\displaystyle=\sum_{j=1}^{N}\frac{a_{j}}{\omega_{j}}I_{j}^{\frac{q}{2}}\bigg{[}\frac{-\cos\tau_{j}}{2m}\bigg{(}\sin^{2m-1}\tau_{j}+\sum_{k=0}^{m-1}s_{k}^{2}\sin^{2m-2k-1}\tau_{j}\bigg{)}+\frac{(2m-1)!!}{2^{m}m!}\tau_{j}\bigg{]}$
$\displaystyle+\sum_{j=1}^{N}\frac{b_{j}}{\omega_{j}}I_{j}^{\frac{q}{2}}\bigg{[}\frac{\sin\tau_{j}}{2m}\bigg{(}\cos^{2m-1}\tau_{j}+\sum_{k=0}^{m-1}s_{k}^{2}\cos^{2m-2k-1}\tau_{j}\bigg{)}+\frac{(2m-1)!!}{2^{m}m!}\tau_{j}\bigg{]},$
(57)

where $s_{k}^{2}=((2m-1)(2m-3)\ldots(2m-2k+1))/(2^{k}(m-1)(m-2)\ldots(m-k))$.
The above expressions can be derived by considering the reduction formulas for
the integrals of $\sin^{n}\theta$ and $\cos^{n}\theta$ and they can be written
as

$\displaystyle\int{\sin^{n}\theta d\theta}$ $\displaystyle=$
$\displaystyle-\frac{1}{n}\sin^{n-1}\theta\cos\theta+\frac{n-1}{n}\int{\sin^{n-2}\theta
d\theta},$ $\displaystyle\int{\cos^{n}\theta d\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\cos^{n-1}\theta\cos\theta+\frac{n-1}{n}\int{\cos^{n-2}\theta
d\theta}.$

##  4 Non-isochronous Systems

In this section, we present the method of generating nonlinear oscillators
having non-isochronous property. Differing from isochronous systems, these
systems have the property of amplitude dependent frequency of oscillations. We
point out below that this property essentially arises from the fact that the
the form of the function $f(x,y)$ should be related to an integral of motion
and we make the simplest of choice $f(x,y)=f(I)$, where $I=x^{2}+y^{2}$, the
energy integral of the linear harmonic oscillator. the details are given in
the following. To make our investigations systematic, we start our discussion
by considering a system of two coupled first order ODEs.

###  4.1 2-Coupled Non-iscochronous Systems

In this case the nonlinear ODEs take the form (see Eq.(9)):

$\displaystyle\dot{x}=\omega f(x,y)y,~{}~{}\dot{y}=-\omega f(x,y)x.$ (58)

From Eq. (58) we find

$\displaystyle\frac{y\dot{x}-x\dot{y}}{(x^{2}+y^{2})}=f(x,y)\omega.$ (59)

Upon introducing the angle variable expression (13) in Eq. (59), the latter
equation becomes

$\displaystyle\dot{\theta}=\Omega(x,y),$ (60)

where $\Omega=\omega f(x,y)$. One may observe that unlike the isochronous
case, here the angular velocity is not a constant but a function of $x$ and
$y$. This constraint introduces a connection between amplitude and the
frequency, as we see below.

Now we analyze the structure of the form $f(x,y)$. Since it appears inside the
phase it should be a constant. So we can consider any form of a constant
associated with the function $f(x,y)$ and consequently the resultant system
becomes a nonisochronous one. Obviously, the simplest choice one can make for
the function $f(x,y)$ is that it can be a function of an integral of motion,
that is $f(x,y)=f(I)$. This considered integral should be independent of $t$.
Since we are studying a second order system the obvious choice is the first
integral ($I=H=x^{2}+y^{2}$). With these arguments we choose
$f(x,y)=f(H(x,y))$.

To obtain the general solution of (58), we follow the same procedure as
discussed in the case of isochronous systems. Considering $H=x^{2}+y^{2}$ and
using the expression (13), we find

$\displaystyle y^{2}=H-y^{2}\tan^{2}(\Omega t+\delta),$ (61)

from which we can fix

$\displaystyle y=\sqrt{H}\cos(\Omega t+\delta).$ (62)

Now substituting (62) in (13), we obtain

$\displaystyle x=\sqrt{H}\sin(\Omega t+\delta).$ (63)

Replacing $\Omega=\omega f(H)$ in Eqs. (62) and (63) we end up with

$\displaystyle x(t)=\sqrt{H}\sin\big{(}f(H)\omega
t+\delta\big{)},~{}~{}~{}y(t)=\sqrt{H}\cos\big{(}f(H)\omega t+\delta\big{)}.$
(64)

The methodology given above clearly demonstrates the interdependency between
amplitude and the frequency of oscillations in a class of nonlinear non-
isochronous systems.

###  4.2 Generalization to $2N$ coupled first order non-isochronous nonlinear
equations

In this section, we generalize the results given in the previous subsection to
the system of $N$ coupled nonlinear oscillators of Li$\acute{e}$nard type.
Here we consider the equation of motion in the form

$\displaystyle\dot{x}_{i}=f(\overline{x},\overline{y})\omega_{i}y_{i},~{}~{}\dot{y}_{i}=-f(\overline{x},\overline{y})\omega_{i}x_{i},\;\;\;\;i=1,2,\cdots,N,$
(65)

where
$f(\overline{x},\overline{y})=f(x_{1},x_{2},...,x_{N},y_{1},y_{2},...,y_{N})$.
To construct the general solution of Eq. (65), we proceed in the following
way.

From Eq. (65) we can derive the following expressions, namely

$\displaystyle
y_{i}\dot{x}_{i}-x_{i}\dot{y}_{i}=\omega_{i}f(\overline{x},\overline{y})(x_{i}^{2}+y_{i}^{2}),\hskip
28.45274pti=1,2,...,N.$ (66)

Let us choose
$f(\overline{x},\overline{y})=f(\overline{H})=f(H_{1},H_{2},\cdots,H_{N})$,
where

$\displaystyle x_{i}^{2}+y_{i}^{2}=H_{i},\qquad i=1,2,\dots$ (67)

Substituting the form $f(\overline{x},\overline{y})=f(\overline{H})$ in (66),
we obtain

$\displaystyle
y_{i}\dot{x}_{i}-x_{i}\dot{y}_{i}=\omega_{i}f(\overline{H})((x_{i}^{2}+y_{i}^{2}),\hskip
28.45274pti=1,2,...,N.$ (68)

Following the steps given in Sec. 5, we derive the solution of Eq. (65) in the
form

$\displaystyle
x_{i}(t)=\sqrt{H_{i}}\sin\Bigg{(}f(\overline{H})\omega_{i}t+\delta_{i}\Bigg{)},\hskip
28.45274pti=1,2,...,N,$ (69) $\displaystyle
y_{i}(t)=\sqrt{H_{i}}\cos\Bigg{(}f(\overline{H})\omega_{i}t+\delta_{i}\Bigg{)},\hskip
28.45274pti=1,2,...,N,$ (70)

where $\delta_{i},~{}i=1,2,3,...,2N$, are integration constants. Out of the
$2N$ integrals, $N$ integrals show explicit time dependency as

$\displaystyle\delta_{i}=\cot^{-1}\bigg{[}\frac{y_{i}}{x_{i}}\bigg{]}-\Omega_{i}t=\cot^{-1}\bigg{[}\frac{y_{i}}{x_{i}}\bigg{]}-f(\overline{H})\omega_{i}t,\hskip
28.45274pti=1,2,...,N$ (71)

and the remaining $N$ integrals exhibit time independency in the form

$\displaystyle H_{i}=x_{i}^{2}+y_{i}^{2},\hskip 28.45274pti=1,2,...,N.$ (72)

These $2N$ functionally independent integrals provide the general solution of
(65).

###  4.3 Examples: Nonisochronous oscillators

####  4.3.1 Example: Mathews-Lakshmanan Oscillator

To demonstrate the theory presented above, we choose
$f(x,y)=[1+\lambda(x^{2}+y^{2})]^{-\frac{1}{2}}$ so that Eq. (58) becomes

$\displaystyle\dot{x}=\frac{\omega
y}{\sqrt{1+\lambda(x^{2}+y^{2})}},~{}~{}~{}\dot{y}=-\frac{\omega
x}{\sqrt{1+\lambda(x^{2}+y^{2})}}.$ (73)

where $\lambda$ and $\omega^{2}$ are arbitrary constants. Rewriting Eq. (73)
as a second order differential equation in $x$, we obtain the celebrated
Mathews-Lakshmanan oscillator equation laksh

$\displaystyle(1+\lambda x^{2})\ddot{x}-\lambda x\dot{x}^{2}+\omega^{2}x=0.$
(74)

The restriction $\lambda=0$ gives harmonic oscillator equation. Eq.(74) admits
the Lagrangian and Hamiltonian structures, respectively, of the form

$\displaystyle L=\frac{\dot{x}^{2}-\omega^{2}x^{2}}{2(1+\lambda
x^{2})}~{}\mathrm{and}~{}H=\frac{1}{2}[p^{2}(1+\lambda
x^{2})+\frac{\omega^{2}x^{2}}{(1+\lambda x^{2})}],$ (75)

where the canonically conjugate momentum $p$ is given by
$p=\frac{\dot{x}}{(1+\lambda x^{2})}$.

During the past few years several studies have been made to explore various
physical and mathematical properties of the nonlinear oscillator equation
(74), see for example Refs.Venk -newone .

The general solution of (73) can be identified from Eq. (64) in the form

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle A\sin\big{(}(1+\lambda
A^{2})^{-1/2}\omega t+\delta\big{)},$ $\displaystyle y(t)$ $\displaystyle=$
$\displaystyle A\cos\big{(}(1+\lambda A^{2})^{-1/2}\omega t+\delta\big{)},$
(76)

where $A$ and $\delta$ are integration constants. Here $A$ is related to the
energy integral by $H=E=\frac{\frac{1}{2}wA^{2}}{1+\lambda A^{2}}$.

In Figure 3 we demonstrate both the harmonic oscillator Eq. (2) and the
nonlinear oscillator Eq. (73) exhibiting periodic oscillations with the same
amplitude for identical initial conditions. However, the nonlinear oscillator
Eq. (73) exhibits amplitude dependent frequency of oscillations. By changing
the initial conditions, the amplitude of the oscillations also gets changed.
Depending on the value of the amplitude $A$, the frequency $\Omega$ of the
nonlinear oscillator (73) also changes. Further, this is also confirmed by the
solution given in (76), where the amplitude $A$ is related with the frequency
$\Omega$ through the expression $\Omega^{2}=\frac{\omega^{2}}{1+\lambda
A^{2}}$.

Figure 3: Time series plot of equations (2) and (73) (red and black)
exhibiting periodic oscillations with same amplitude but with different
frequencies for different initial conditions (a) $x(0)=u(0)=-0.5$ and
$y(0)=v(0)=0.5$ and (c) $x(0)=u(0)=-1$ and $y(0)=v(0)=0.5$. Corresponding
phase space portraits are given in (b) and (d) for $\lambda=\omega=1.0$.

####  4.3.2 Example: 2-coupled nonlinear Mathews-Lakshmanan (ML) oscillator
equation

Here, we consider a system of 2-coupled nonlinear ML oscillator equation whose
equation of motion is given by

$\displaystyle\dot{x}_{1}=fw_{1}y_{1},~{}~{}\dot{x}_{2}=fw_{2}y_{2},~{}~{}\dot{y}_{1}=-fw_{1}x_{1},~{}~{}\dot{y}_{2}=-fw_{2}x_{2},$
(77)

where $f=[1+\lambda(x_{1}^{2}+x_{2}^{2}+y_{1}^{2}+y_{2}^{2})]^{-\frac{1}{2}}$.
Rewriting the above four first order ODEs into two coupled second order ODEs
in the variables $x_{1}$ and $x_{2}$, we find

$\displaystyle(1+\lambda(x_{1}^{2}+x_{2}^{2}))\ddot{x}_{1}-\lambda
x_{1}(\dot{x}_{1}^{2}+\frac{w_{1}^{2}}{w_{2}^{2}}\dot{x}_{2}^{2})+w_{1}^{2}x_{1}=0,$
$\displaystyle(1+\lambda(x_{1}^{2}+x_{2}^{2}))\ddot{x}_{2}-\lambda
x_{2}(\frac{w_{2}^{2}}{w_{1}^{2}}\dot{x}_{1}^{2}+\dot{x}_{2}^{2})+w_{2}^{2}x_{2}=0,$
(78)

where $\lambda$, $w_{1}$ and $w_{2}$ are arbitrary parameters.

Following the procedure given in the previous section, we find the general
solution in the form

$\displaystyle
x_{1}=\sqrt{H_{1}}\sin((1+\lambda(H_{1}+H_{2})^{(-1/2)})w_{1}t+\delta_{1}),$
$\displaystyle
y_{1}=\sqrt{H_{1}}\cos((1+\lambda(H_{1}+H_{2})^{(-1/2)})w_{1}t+\delta_{1}),$
$\displaystyle
x_{2}=\sqrt{H_{2}}\sin((1+\lambda(H_{1}+H_{2})^{(-1/2)})w_{2}t+\delta_{2}),$
$\displaystyle
y_{2}=\sqrt{H_{2}}\cos((1+\lambda(H_{1}+H_{2})^{(-1/2)})w_{2}t+\delta_{2}),$
(79)

where $H_{i}$ and $\delta_{i},~{}i=1,2,3,4$, are integration constants. One
important point we wish to note here is that the general solution for the
coupled ML oscillators has been given for the case $w_{1}\neq w_{2}$. In Ref.
sir1 , the authors have generalized the one dimensional ML oscillator to two
dimensions as well as $n$-dimensions by considering the parameter $w$ as same
in all the dimensions.

####  4.3.3 Example: $2N$-coupled nonlinear ML oscillator equation

The equation of motion for the $2N$-coupled nonlinear ML oscillator can be
written as

$\displaystyle\dot{x}_{i}=f(\overline{x},\overline{y})\omega_{i}y_{i},~{}~{}\dot{y}_{i}=-f(\overline{x},\overline{y})\omega_{i}x_{i},\;\;\;\;i=1,2,\cdots,N,$
(80)

where
$f(\overline{x},\overline{y})=f=[1+\lambda(x_{1}^{2}+x_{2}^{2}+\cdots+x_{N}^{2}+y_{1}^{2}+y_{2}^{2}+\cdots+y_{N}^{2})]^{-\frac{1}{2}}$.
Rewriting the above $2N$ first order ODEs into $N$-coupled second order ODEs
in the variables $x_{1},x_{2},\cdots,x_{N}$, we find

$\displaystyle(1+\lambda\sum_{j=1}^{N}x_{j}^{2})\ddot{x}_{i}-\lambda
x_{i}\sum_{j=1}^{N}\frac{w_{i}^{2}}{w_{j}^{2}}\dot{x}_{j}^{2}+\omega_{i}^{2}x_{i}=0,$
(81)

where $\lambda$ and $\omega_{i}^{2}$ are arbitrary constants. Choosing
$f(\overline{x},\overline{y})=f(\overline{H})=(1+\lambda\sum_{i=1}^{N}H_{i})^{-1/2}$,
we obtain the general solution of $(\ref{cmod96})$ in the form

$\displaystyle x_{i}(t)$ $\displaystyle=$
$\displaystyle\sqrt{H_{i}}\sin\Bigg{(}\bigg{(}1+\lambda\sum_{j=1}^{N}H_{j}\bigg{)}^{-1/2}\omega_{i}t+\delta_{i}\Bigg{)},$
$\displaystyle y_{i}(t)$ $\displaystyle=$
$\displaystyle\sqrt{H_{i}}\cos\Bigg{(}\bigg{(}1+\lambda\sum_{j=1}^{N}H_{j}\bigg{)}^{-1/2}\omega_{i}t+\delta_{i}\Bigg{)},\hskip
2.84544pti=1,2,...,N,$ (82)

where $H_{i}$ and $\delta_{i}$ are the integration constants.

##  5 Conclusion

In this work, we have discussed a method of generating isochronous and non-
isochronous nonlinear oscillators from the integrals of the simple harmonic
oscillator equation. By considering the integrals to be in the same form for
both the linear and nonlinear oscillators, we have identified the nonlinear
oscillators that posses either isochronous and non-isochronous property
(involving harmonic functions). We then considered two coupled first order
ODEs and discussed the method of constructing the solution for both the cases.
In each case we have also demonstrated the theory with examples. We then
extended the theory to $2N$ coupled first order ODEs. We have also derived the
general solution and presented the explicit forms of integrals for both the
cases. We have also given examples for each of the cases. As far as the
isochronous systems are concerned, the value of the exponent $q$ for the
function $f$ in Sec. 3 decides the nature of the solution which will be either
periodic or damped oscillations. The general solution and the integrals for
the $2N$-coupled ML oscillator system which possesses the non-isochronous
property are also presented for the first time in the literature.

##  6 Acknowledgments

RMS is funded by the Center for Nonlinear Systems, Chennai Institute of
Technology, India, vide funding number CIT/CNS/2021/RP-015. V.K.C. thanks DST,
New Delhi for computational facilities under the DST-FIST programme (Grant No.
SR/FST/PS- 1/2020/135) to the Department of Physics. V.K.C. is also supported
by the SERB-DST-MATRICS Grant No. MTR/2018/000676. The work of M.S forms part
of a research project sponsored by National Board for Higher Mathematics
(NBHM), Government of India under Grant No. 02011/20/2018NBHM(R.P)/R&D
II/15064 and M.L. is supported by a DST-SERB National Science Chair.

## References

  * (1) R. Iacono and F. Russo, Class of solvable nonlinear oscillators with isochronous orbits. Phys. Rev. E 83, (2011) 027601.
  * (2) V. K. Chandrasekar, Jane H. Sheeba, R. Gladwin Pradeep, R. S. Divyasree and M. Lakshmanan, A class of solvable coupled nonlinear oscillators with amplitude independent frequencies. Phys. Letter A 376, (2012) 2188-2194.
  * (3) U. Tanriver, S. Roy Choudhury and G. Gambino, Lagrangian dynamics and possible isochronous behavior in several classes of nonlinear second order oscillators via the use of Jacobi last multiplier. Int. J. Nonlinear Mech. 74, (2015) 100-107.
  * (4) A. G. Choudhury and P. Guha, On isochronous cases of the Cherkas system and Jacobi’s last multiplier. J. Phys. A: Math. Theor. 43, (2010) 125202.
  * (5) J. F. Cari$\tilde{n}$ena, A. M. Perelomov and M. F. Ra$\tilde{n}$ada, Isochronous classical systems and quantum systems with equally spaced spectra. J. Phys.: Conf. Ser. 87 (2007) 012007.
  * (6) F. Calogero, Isochronous dynamical systems. Phil. Trans. R. Soc. A 369, (2011) 1118-1136.
  * (7) V. M. Budanov, On a nonlinear isochronous system. Moscow Univ. Mech. Bull. 68, (2013) 145.
  * (8) A. G. Choudhury, P. Guha. Chiellini integrability condition, planar isochronous systems and Hamiltonian structures of Li$\acute{e}$nard equation. Discrete and Continuous Dynamical Systems - B, 22 (2017) 2465-2478.
  * (9) P. Guha and A. G. Choudhury, The Jacobi last multiplier and isochronicity of Li$\acute{e}$nard type systems. Rev. Math. Phys., 25 (2013), 1330009.
  * (10) A. Algaba, E. Freire and E. Gamero, Isochronicity via normal form. Qual. Theory Dyn. Syst. 1, 2 (2000), 133-156.
  * (11) M. Bardet, I. Boussaada, A. R. Chouikha and J. M. Strelcyn, Isochronicity conditions for some planar polynomial systems II. Bull. Sci. math. 135 (2011) 230-249.
  * (12) O. A. Chalykh and A. P. Veselov, A remark on rational isochronous potentials. J. Nonlinear Math. Phys. 12 (2005) 179-183.
  * (13) R. Cobiaga and W. Reartes, Isochronous families of limit cycles. Electronic J. Differ. Eq., 117, (2018), 1-16.
  * (14) V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, A unification in the theory of linearization of second-order nonlinear ordinary differential equations, J. Phys. A 39 (2006) L69-L76.
  * (15) L. G. S. Duarte, C. Moreira and F. C. Santos, Linearization under nonpoint transformations, J. Phys. A 27 (1994) L739-L743.
  * (16) C. Muriel and J. L. Romero, Nonlocal transformations and linearization of second-order ordinary differential equations, J. Phys. A: Math. Theor. 43 (2010) 434025.
  * (17) A. Durga Devi, R. Gladwin Pradeep, V. K. Chandrasekar and M. Lakshmanan, Method of generating N-dimensional isochronous nonsingular Hamiltonian systems. J. Nonlinear Math. Phys. 104, (2017) 78-93.
  * (18) J. Chavarriga and M. Sabatini, A survey of isochronous centers. Qual. Theory Dyn. Syst. 1, (1999) 1-70.
  * (19) A. R. Chouikha, Isochronous centers of Li$\acute{e}$nard type equations and applications. J. Math. Anal. Appl., 331 (2007), 358-376.
  * (20) F. Calogero, Isochronous systems. Oxford University Press, Oxford, (2008).
  * (21) F. Calogero and F. Leyvraz, How to extend any dynamical system so that it becomes isochronous, asymptotically isochronous or multi-periodic. J. Nonlin. Math. Phys. 16, (2009) 311-338.
  * (22) F. Calogero and F. Leyvraz, General technique to produce isochronous Hamiltonians. J. Phys. A 40, (2007) 931-944.
  * (23) F. Calogero and F. Leyvraz, A new class of isochronous dynamical systems. J. Phys. A: Math. Theor. 41 (2008) 295101.
  * (24) P. Guha and A. G. Choudhuri, Symplectic rectification and isochronous Hamiltonian systems. J. Phys. A 42, (2009) 192001.
  * (25) M. Lakshmanan and S. Rajasekar, Nonlinear dynamics-Integrability, chaos and patterns, Springer-Verlag, Berlin Heidelberg, (2003).
  * (26) P. M. Mathews and M. Lakshmanan, On a unique nonlinear oscillator, Quart. Appl. Math. 32, (1974) 215-218.
  * (27) P. M. Mathews and M. Lakshmanan, A quantum mechanically solvable nonpolynomial Lagrangian with velocity-dependent interaction, Nuovo Cimento A 26, (1975) 299-315.
  * (28) H. Daido and K. Nakanishi, Diffusion-induced inhomogeneity in globally coupled oscillators: Swing-by mechanism, Phys. Rev. Lett. 96, (2006) 054101.
  * (29) M. Sebek, R. T$\ddot{o}$njes, and I. Z. Kiss, Complex rotating waves and long transients in a ring network of electrochemical oscillators with sparse random cross-connections, Phys. Rev. Lett. 116, (2016) 068701.
  * (30) M. Wickramasinghe and I. Z. Kiss, Synchronization of electrochemical oscillators with differential coupling, Phys. Rev. E 88, (2013) 062911.
  * (31) I. S. Gradshteyn and I. M. Ryzhik, Table of integrals, series, and products, (Editors: A. Jeffrey and D. Zwillinger), Academic Press, San Diego (2000).
  * (32) V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, Unusual Li$\acute{e}$nard-type nonlinear oscillator, Phy. Rev. E 72, (2005) 066203.
  * (33) W. S. Loud, Behavior of the periodic of solutions of certain plane autonomous systems near centers, Contributions to Differential Equations. 3, (1964) 21-36.
  * (34) S. Karthiga, V. Chithiika Ruby, M. Senthilvelan, and M. Lakshmanan, Quantum solvability of a general ordered position dependent mass system: Mathews-Lakshmanan oscillator, J. Math. Phys. 58 (2017) 102110.
  * (35) V. Chithiika Ruby and M. Lakshmanan, Quantum solvability of quadratic Li$\acute{e}$nard type nonlinear oscillators possessing maximal Lie point symmetries: An implication of arbitrariness of ordering parameters, J. Phys. Commun. 5 (2021) 065007.
  * (36) A. Bhuvaneswari, V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, On the complete integrability of a nonlinear oscillator from group theoretical perspective, J. Math. Phys. 53, (2012) 073504.
  * (37) J. F. Cari$\tilde{n}$ena, J. de Lucas J, M. F. Ra$\tilde{n}$ada, Jacobi multipliers, non-local symmetries and nonlinear oscillators, J. Math. Phys. 56, (2015) 063505.
  * (38) O. Mustafa, $n$-dimensional PDM-damped harmonic oscillators: linearizability, and exact solvability, Phys. Scr. 96, (2021) 065205.
  * (39) O. Mustafa, Isochronous $n$-dimensional nonlinear PDM-oscillators: linearizability, invariance and exact solvability. Eur. Phys. J. Plus 136, (2021) 249.
  * (40) M. Lakshmanan and V. K. Chandrasekar, Generating finite dimensional integrable nonlinear dynamical systems, Eur. Phys. J. Special Topics. 222, (2013) 665-688.
  * (41) J. F. Cari$\tilde{n}$ena, M. F. Ra$\tilde{n}$ada, M. Santander and M. Senthilvelan, A non-linear oscillator with quasi-harmonic behaviour: two- and $n$-dimensional oscillators, Nonlinearity. 17, (2004) 1941-1963.