# Tailoring radiation pressure on infinite slab using pair of non-collinear
plane waves

R. Ali<EMAIL_ADDRESS>Applied Physics Department, Gleb Wataghin
Physics Institute, University of Campinas, Campinas 13083-859, SP, Brazil  R.
S. Dutra  LISComp-IFRJ , Instituto Federal de Educação, Ciência e Tecnologia,
Rua Sebastião de Lacerda, Paracambi, RJ, 26600-000, Brasil  S. Iqbal
Department of Physics, School of Natural Sciences National University of
Sciences and Technology Sector H-12, Islamabad, Pakistan

###### Abstract

The electromagnetic field exerts radiation pressure on the matter and tends to
move it either in the backward or forward direction due to net optical pulling
or pushing force, respectively. In this work, we reveal an interesting
phenomenon of a local positive and negative radiation pressure on a dielectric
(chiral) slab by using two linearly (circularly) polarized plane waves. In
this regard, we develop for the first time, a theory to describes the local
radiation pressure appearing due to the interference between the two obliquely
impinging (non-collinear) light sources. Under this situation, the radiation
pressure strongly depends on the angle of incidence, the polarization of the
electromagnetic field and the chirality parameters of the slab (in the case of
chiral medium). Our numerical analysis shows that the radiation pressure,
exerted on a dielectric or a chiral slab due to the two incident plane waves,
are constant over the slab for normal incidence, and it varies locally for an
oblique incidence, which indeed follows the conservation laws at all incident
angles. It is expected that the results may find fruitful applications in
optical manipulation of soft matters, for instance, cell membranes, chiral
surfaces and other soft materials.

##  I Introduction

The electromagnetic field carries both linear and angular momentum which can
be exchanged during the field-matter interaction, and yield a radiation
pressure and radiation torque on the matter Ashkin1970 . In the past few
decades, the field of light-matter interaction has grown impressively after
the invention of optical tweezers ashkin1986  that have played an
indispensable role with countless applications Gieseler2020 ; Polimeno2018 .
In optical tweezers setups, the trapping and manipulating of spherical objects
are obtained near to the focus point of a tightly focused continuous-wave
laser beam ashkin2006 ; Friese1998 ; grier2003 ; Lin2014 , while the long-
range optical manipulation of chiral and dielectric objects can be achieved by
using a structured Bessel beam or two non-collinear plane waves Ali2020Tailor
; Ali2021 ; chan2011 ; li2019 ; Shvedov2014 ; dogariu2011 .

The theory of optical manipulation is not limited to spherical particles but
tremendous efforts have been made to calculate the radiation pressure on
various kinds of objects, such as, composite spheres Ali2020Tailor ; Ali2018 ,
spheroidal particles Stephen2007 , nanorices Wiley2007 , cylindrical particles
Gauthier1999 , chiral spheres Ali2020 ; Ali2020josa ; Canaguier2013 ; Kun2014
, chiral nanocrystals Ali2020  and dielectric slabs Mansuripur2005 ;
Loudon2010 ; Alaee2018 , wherein the magnitude as well as direction of the
optical forces and the flow of optical momenta have been discussed. Nowadays,
metamaterials and metasurfaces are emerging fields that have potential
applications in the field of nanotechnology, sensing and medicine
Corradini2007 . The negative radiation pressure on such surfaces, for
instance, cell membranes, and chiral metasurfaces Basiri  could be used to
probe the local properties by stretching and relaxing the surfaces due to the
incident light. To the best of our knowledge, the local radiation pressure on
such metasurfaces has not been discussed so far.

It is worth mentioning that a single plane wave is obliged to exert positive
radiation pressure on any passive system in the paraxial approximation due to
momentum conservation Lepeshov2020 ; Ali2020 ; ashkin2006 . For example, the
optical force $F$ acting on a slab located in vacuum can be expressed as
$F=N_{i}\hbar k{(1+r-t)}$, where
$N_{i}=c\varepsilon_{0}E_{0}^{2}/2\hbar\omega$ is the photon density,
$E_{0}^{2}$ is the electric field intensity, c is the speed of light, $\omega$
is the angular frequency of the light and $\varepsilon_{0}$ is the vacuum
permittivity. Furthermore, $t$ and $r$ are the Fresnel transmission and
reflection coefficients that use to define the momentum transmitted and
reflected through the slab, respectively. In this regard, the optical force
exerted on a perfect reflector  ($r=1,\,t=0$), black body (${r}=0,\,t=0$), and
perfect transmitter ($r=0,\,t=1$) can be estimated as $F=2N_{i}\hbar k$,
$F=N_{i}\hbar k$ and $F=0$, respectively. Thus, the net force can not be
negative due to fact that the transmitted or reflected photons cannot have
larger momentum than the incident ones and hence, $F\geq 0$ Lepeshov2020 ;
Mizrahi2010 ; Wang2015 ; Li2016 .

However, the situation changes dramatically when we employ optical gain
Ali2021gain  or shine (normally) a structured incident beam such as Bessel
beam. In such cases, one can tune the direction of scattered momentum flow in
forward direction and due to conservation of momentum, the recoiling force
could be negative. In all the cases, the optical force or indeed the radiation
pressure is constant over the slab, thereby, net force appears at the center
of mass of the slab.

In this paper, we are intended to calculate the locally negative and positive
radiation pressure density exerted by a pair of plane waves on the slab by
evaluating the Maxwell stress tensor. To this end, we consider an oblique
incidence of two linearly polarized plane waves separated by an angular
distance $2\theta$ (see Sec. II ) impinging on the dielectric slab.
Furthermore, we extend this framework to calculate the negative radiation
pressure exerted by two circularly polarized non-collinear plane waves on
chiral slabs, for instance, slabs wherein the chirality may appear due to
randomness, chiral metasurfaces Basiri  or naturally occurring biological
surfaces. In all the cases, we examine the radiation pressure acting on the
slab due to interference of incident plane waves and analyze the role of
incident angle $\theta$ and chirality parameter. In this perspective, we
presented for the first, a new analytical framework to calculate the momentum
density over the slab due to the superposition of plane waves of light, which
meets all the necessary conservation laws for all the incident angles.

The rest of the paper is organized as follows. Section II is devoted to the
derivation of our theoretical formalism, where we derive the optical radiation
pressure expressions for dielectric and chiral slab and discussed its limiting
cases to well-known results. In Sec. III, we present numerical results and
analyze the findings. Finally, we summarize our findings and conclusions in
Sec IV.

##  II Methodology

###  II.1 Radiation pressure on dielectric slab due to two linearly polarized
plane waves

We consider two linearly polarized transverse electric (TE) plane waves
impinging obliquely at an angle $\theta_{i}$ (where $i=1,2$) on an infinite
dielectric slab of thickness $d$ and refractive index
$n_{s}=\sqrt{\varepsilon_{s}\mu_{s}}$, where $\varepsilon_{s}$ and $\mu_{s}$
are the relative permittivity and relative permeability of the slab, as
depicted in figure 1. The total impinging electric field on the slab is
defined by taking the superposition of both plane waves propagating through
the host medium of relative permittivity $\varepsilon$ as

$\mathbf{E}_{in,total}=E_{0}\sum_{j=1}^{2}e^{i\mathbf{k}_{j}\cdot\mathbf{r}}\hat{y},$
(1)

where $\mathbf{k}_{j}$ ($j=1,2$) is the wave vector of $j^{th}$ plane wave and
can be expressed as
$\mathbf{k}_{j}=(-1)^{j+1}k\sin\theta_{j}\mathbf{\hat{x}}+k\cos\theta_{j}\mathbf{\hat{z}}$.
At the boundary of the slab, a part of $\mathbf{E}_{in,total}$ will be
reflected back into surrounding medium and propagate with wave vector
$\mathbf{k}^{\prime}_{j}=(-1)^{j+1}k\sin\theta\mathbf{\hat{x}}-k\cos\theta\mathbf{\hat{z}}$,
for the sake of easiness we take $\theta_{1}=\theta_{2}=\theta$. By conserving
the polarization, the total reflected field $\mathbf{E}_{r,total}$ can be
written as

Figure 1: Schematic diagram: an infinite dielectric slab of thickness $d$,
relative permittivity $\epsilon_{s}$ is illuminated by two non-collinear
transverse electric plane waves making angle $\theta$ with the z-axis, where
dashed lines show the surfaces of the imaginary Gaussian surfaces with normal
unitary vector $\hat{n}.$ Radiation pressure acting on the slab interface
$z=0$ (green) and $z=d$ (dotted) and total local radiation pressure on
$F({x,\theta})=F({z=0,x,\theta})+F({z=d,x,\theta})$ (black) for (b) normal
incidence and (c) oblique incidence at $\theta=70^{o}$.

$\mathbf{E}_{r,total}=E_{0}\,r_{TE}\sum_{j=1}^{2}e^{i\mathbf{k}^{\prime}_{j}\cdot\mathbf{r}}\hat{y},$
(2)

where $r_{TE}$ is reflection coefficient for infinite slab Born . The rest
part of $\mathbf{E}_{in,total}$ is transmitted through the slab at $z=d$ with
the same incident angle $\theta$ and allowing us to write the total
transmitted electric field as

$\mathbf{E}_{t,total}=E_{0}\,t_{TE}\sum_{j=1}^{2}e^{i\mathbf{k}_{j}\cdot\mathbf{r}}\hat{y},$
(3)

where, $t_{TE}$ is the transmission coefficient of an infinity slab Born .

Likewise, the similar expressions for the total incident, reflected and
transmitted magnetic fields can also be obtained by using the relation
$\mathbf{H}_{j}=\frac{\mathbf{k}_{j}}{k_{j}\mu v}\times\mathbf{E}_{j}$ for
each plane wave, where $\mu$ and $v$ are the host medium magnetic permeability
and the speed of light, respectively.

The radiation pressure on an infinite slab is calculated by using the Maxwell
stress tensor $\overleftrightarrow{T}$ Brandon  through the closed Gaussian
surfaces as represented by the dashed lines in Fig. 1 Brandon  on the both
sides of the slab, such as $z=0^{-}$ and $z=d^{+}$, with respective normal
unitary vectors $\mathbf{\hat{n}}=\mathbf{-\hat{z}}$ and
$\mathbf{\hat{n}}=\mathbf{\hat{z}}$, respectively. Finally, the radiation
pressure $\mathbf{F}$ can be defined as

$\mathbf{F}=\frac{1}{2}\Re\\{\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=0^{-})-\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=d^{+})\\},$
(4)

where $\overleftrightarrow{T}(z=0)$ and $\overleftrightarrow{T}(z=d^{+})$ are
the Maxwell stress tensor evaluated at point z =0 and z=d. Without going to
the debate on the Maxwell stress tensors, here we are using Minkowski Maxwell
stress tensor that take into account the medium contribution Pfeifer2007 . The
explicit expression for $\overleftrightarrow{T}$ can be found by substituting
the total electric and magnetic fields in the standard expression (see Eq. 4
of the Ref. Brandon ) of the Maxwell stress tensor and is expressed as

$\displaystyle\overleftrightarrow{T}=\frac{I}{2}(\varepsilon\mathbf{E}_{total}\cdot\mathbf{E}^{*}_{total}+\mu\mathbf{H}_{total}\cdot\mathbf{H}^{*}_{total})$
$\displaystyle-\varepsilon\mathbf{E}_{total}\mathbf{E}^{*}_{total}-\mu\mathbf{H}^{*}_{total}\mathbf{H}_{total}$
(5)

where
$\mathbf{E}_{total}(z=0^{-})=\mathbf{E}_{in,total}(z=0^{-})+\mathbf{E}_{r,total}(z=0^{-})$,
$\mathbf{E}_{total}(z=d^{+})=\mathbf{E}_{t,total}(z=d^{+})$ and $\varepsilon$
is relative permittivity of the surrounding medium. By setting the Eqs.
(1)-(3), we obtain the total radiation pressure on the slab by two non-
collinear plane waves

$\displaystyle\mathbf{F}=\mathbf{F}(x,\theta)=\varepsilon
E_{0}^{2}\biggl{[}\biggl{(}|1+r_{TE}|^{2}+|1-r_{TE}|^{2}\cos^{2}\theta$
$\displaystyle-|t_{TE}|^{2}(1+cos^{2}\theta)\biggr{)}\,\cos^{2}(kx\sin\theta){+}\biggl{(}|t_{TE}|^{2}-|1+r_{TE}|^{2}\biggl{)}$
$\displaystyle\sin^{2}\theta\,\sin^{2}(kx\sin\theta)\biggr{]}\mathbf{\hat{z}}.$
(6)

Since the problem involves the electromagnetic field solutions inside and
outside the slab and the transmission $t_{TE}$ and reflection $r_{TE}$
coefficients are obtained by applying the appropriate boundary conditions at
the interfaces between the surrounding medium and the slab. The final
expression are given as

$r_{TE}=\frac{r_{m,s}+r_{s,m}exp(i2\beta)}{1+r_{m,s}r_{s,m}exp(i2\beta)}$ (7)

and

$t_{TE}=\frac{t_{m,s}t_{s,m}exp(i\beta)}{1+r_{m,s}r_{s,m}exp(i2\beta)},$ (8)

where, $\beta=2\pi n_{s}d\cos\theta/\lambda$ is optical path inside the slab.
The coefficients ($t_{m,s}$, $r_{m,s}$) and ($t_{s,m}$, $r_{s,m}$) are
evaluated at the interfaces and are expressed as

$t_{m,s}=\frac{2n_{m}\cos\theta_{m}}{n_{m}\cos\theta_{m}+n_{s}\cos\theta_{s}},$
(9)

$t_{s,m}=\frac{2n_{s}\cos\theta_{s}}{n_{s}\cos\theta_{s}+n_{m}\cos\theta_{m}}$
(10)

and

$r_{m,s}=-r_{s,m}=\frac{n_{m}\cos\theta_{m}-n_{s}\cos\theta_{s}}{n_{m}\cos\theta_{m}+n_{s}\cos\theta_{s}},$
(11)

The refraction angle $\theta_{s}$ inside slab and the incident angle
$\theta_{m}=\theta$ are connected by Snell’s law
$\theta_{s}=\arcsin({n_{m}\sin\theta/n_{s}})$.

In order to gain a more physical insight, we compare the limiting case of our
exact analytical expressions presented in Eq. 6, when $\theta=0$ we recover
the well known results for the radiation pressure acting on an infinite slab
with gain/loss due to a normal incidence Alaee2018 ; Lepeshov2020 ;
Mizrahi2010 ; Brandon .

$\mathbf{F}{{(0\degree)}}=2\varepsilon
E_{0}^{2}\biggl{[}1+|r_{TE}|^{2}-|t_{TE}|^{2}\biggl{]}\mathbf{\hat{z}}.$ (12)

It is clear from Eq. 12, the radiation pressure due to normal incidence does
not depend on transversal position. Thus, the pressure obliged to be constant
on each interface of the slab Fig.1(b). However, when both plane waves are
impinging obliquely then the reflected and transmitted light are expected to
produce an interference pattern along x-direction. Thus, each surface of the
slab undergoes to the local maximum and minimum radiation pressure as shown in
Fig.1(c). In Fig.1(b)-(c) we calculate the radiation pressure on each
interface for fixed thickness $d=\lambda$, $n_{m}=1.332$ and $n_{s}=1.5$. The
details of analytical calculation of the radiation pressure acting on
interfaces $z=0:$ $\mathbf{F}(z=0,x,\theta)$ and $z=d:$
$\mathbf{F}(z=d,x,\theta)$ are presented in appendix A.

###  II.2 Radiation pressure on chiral slab due to two circularly polarized
plane waves

Propagation of the electromagnetic wave through a chiral media can be
explained by the following constitutive relations that connect the complex
displacement field D and magnetic field B Ali2020 ; Ali2020josa ; wang2014 ;
Bohren ; Ali2020Theory  as

$\displaystyle\begin{bmatrix}{\bf D}\\\ {\bf
B}\end{bmatrix}=\begin{bmatrix}\varepsilon_{0}\varepsilon&i\sqrt{\varepsilon_{0}\mu_{0}}\,\kappa\\\
-i\sqrt{\varepsilon_{0}\mu_{0}}\,\kappa&\mu_{0}\mu\end{bmatrix}\begin{bmatrix}{\bf
E}\\\ {\bf H}\end{bmatrix},$ (13)

through a direct coupling between the electric field ${\bf E}$ and the
auxiliary field ${\bf H}$ by chirality parameter $\kappa$. Here, $\varepsilon$
and $\mu$ are the relative permittivity and relative permeability of the
medium, respectively. The constant $\kappa$ characterizes the nature of media
and usually satisfies the condition $\kappa\ll\sqrt{\varepsilon\mu}$. By using
the aforementioned constitutive relations in the Maxwell equations Bohren ,
one can derive the insightful information such that the each circular
polarization field propagates with different phase velocity and wave vectors
that explicitly depends on chiral parameter $\kappa$. The wave vectors inside
the chiral media can be written as
$k_{s,\sigma}=k_{0}(\sqrt{\epsilon_{s}}+\sigma\kappa)$, where $k_{0}=2\pi
c/\lambda_{0}$ the vacuum wave vector, $\sigma=1(-1)$ for left (right)
circularly polarized plane wave.

In this section, we consider two oblique circularly polarized incident plane
waves impinging on an infinite chiral slab, immersed in water as depicted in
Fig. 2. Following the same procedure, as adopted in the previous section, one
can write the total incident field as

$\mathbf{E}_{in,total}=E_{0}\sum_{j=1}^{2}e^{i\mathbf{k}_{j}\cdot\mathbf{r}}\frac{(\mathbf{\hat{p}}_{j}+i\sigma_{j}\mathbf{\hat{s}}_{j})}{\sqrt{2}},$
(14)

with unitary vectors
$\mathbf{\hat{p}}_{j}=cos\,\theta_{j}\,\mathbf{\hat{x}}+(-1)^{j}\,sin\,\theta_{j}\,\mathbf{\hat{z}}$
and $\mathbf{\hat{s}}_{j}=\mathbf{\hat{y}}$. The plane waves impinge the slab
surface with wave vectors
$\mathbf{k}_{j}=(-1)^{j+1}k\sin\theta_{j}\mathbf{\hat{x}}+k\cos\theta_{j}\mathbf{\hat{z}}$,
and we assume, $\theta_{1}=\theta_{2}=\theta$, for simplicity.

The reflected plane waves that propagate with the wave vectors
$\mathbf{k}^{\prime}_{j}=(-1)^{j+1}k\sin\theta_{j}\mathbf{\hat{x}}-k\cos\theta_{j}\mathbf{\hat{z}}$
are expanded in terms of the linearly polarized basis as

Figure 2: Schematic illustration of the system, an infinite slab of thickness
$d$, relative permittivity $\epsilon_{s}$, chirality parameter $\kappa$
illuminated by two non-collinear circularly polarized plane waves impinging at
an angle $\theta$ with respect to the z axis, $\hat{n}$ is the normal unitary
vector.

$\mathbf{E}_{r,total}=\sum_{j=1}^{2}e^{i\mathbf{k}^{\prime}_{j}\cdot\mathbf{r}}\frac{(E_{r,\parallel}\mathbf{\hat{p}}^{\prime}_{j}+E_{r,\perp}\mathbf{\hat{s}}^{\prime}_{j})}{\sqrt{2}},$
(15)

with the reflected unit vectors
$\mathbf{\hat{p}}^{\prime}_{j}=-\cos\theta_{j}\,\mathbf{\hat{x}}+(-1)^{j}\,\sin\theta_{j}\,\mathbf{\hat{z}}$
and $\mathbf{\hat{s}}^{\prime}_{j}=\mathbf{\hat{s}}_{j}=\mathbf{\hat{y}}$.

After crossing the slab, the plane waves are transmitted to dielectric host
medium, with the same incident angles $\theta_{1}=\theta_{2}=\theta$. We also
expand the transmitted total electric field on the linearly polarized basis as

$\mathbf{E}_{t,total}=\sum_{j=1}^{2}e^{i\mathbf{k}_{j}\cdot\mathbf{r}}\frac{(E_{t,\parallel}\mathbf{\hat{p}}_{j}+E_{t,\perp}\mathbf{\hat{s}}_{j})}{\sqrt{2}}.$
(16)

The analogous expressions for the total incident, reflected and transmitted
magnetic fields can be obtained in a trivial way by using the relation
$\mathbf{H}_{j}=\frac{\mathbf{k}_{j}}{k_{j}\mu v}\times\mathbf{E}_{j}$.

Figure 3: Normalized optical radiation pressure acting on a dielectric slab by
plane wave shining at a normal incidence (i. e. we take $\theta=0\degree$) as
a function of $d/\lambda$. The refractive index of the slab is (a)
$\rm{Re}(n_{s})=3.54$ and (b) $\rm{Re}(n_{s})=1.5$ for different imaginary
parts: such as $\rm{Im}(n_{s})=-0.05$ (black), ${\rm Im}(n_{s})=0$ (dashed),
and ${\rm Im}(n_{s})=0.05$ (red).

The next step is to find the reflections ($E_{r,\parallel},E_{r,\perp}$) and
transmissions ($E_{t,\parallel},E_{t,\perp}$) coefficients outside the chiral
slab. For this purpose, the incident electric and magnetic fields defined in
Eqs. (15)-(16), together with expressions for the electromagnetic field inside
the chiral slab, obtained by bassiri , can further be written in different
regions by applying the appropriate boundary conditions. For instance, one can
define: (i) the reflected and transmitted fields inside the host medium and
inside the chiral slab, respectively, at $z=0$, (ii) the reflected and
transmitted fields inside the slab and inside the surrounding medium at $z=d$,
respectively. This is simple but lengthy process and the all calculations are
rigorously presented in the Ref. bassiri . By following the same procedure we
come up with eight equations and eight variables.

For our purpose we solve the reflections ($E_{r,\parallel},E_{r,\perp}$) and
transmissions ($E_{t,\parallel},E_{t,\perp}$) coefficients by solving the
system of equations, in the form of $8\times 8$ Matrix, (by following the
procedure as given in Ref. bassiri ) by using the Mathematica. Furthermore, we
consider the both incident plane waves with the same helicities
$\sigma_{1}=\sigma_{2}=\sigma=\pm 1$

The radiation pressure on an infinite slab is determined by solving the
Maxwell stress tensor through closed Gaussian surfaces (represented by the
dashed lines in Fig. 2) Brandon  and given as

$\displaystyle\mathbf{F}=\frac{\varepsilon
E_{0}^{2}}{2}\biggl{[}\biggl{(}|1-E_{r,\parallel}|^{2}+|1+i\sigma
E_{r,\perp}|^{2}-|E_{t,\parallel}|^{2}-|E_{t,\perp}|^{2}\bigg{)}$
$\displaystyle
cos^{2}\theta\,cos^{2}(kx\sin\theta)+\biggl{(}|E_{t,\parallel}|^{2}+|E_{t,\perp}|^{2}-|1+E_{r,\parallel}|^{2}-|1$
$\displaystyle-i\sigma
E_{r,\perp}|^{2}\biggr{)}\biggl{\\{}\sin^{2}\theta\,\sin^{2}(kx\sin\theta)-\cos^{2}(kx\sin\theta)\biggr{\\}}\biggr{]}\mathbf{\hat{z}}.$
(17)

where $\mathbf{F}=\mathbf{F}(x,\theta)$. In this case, the expressions for the
transmission and reflection coefficients are lengthy and complicated, we solve
them by using the Mathematica software. The Mathematica file containing the
solution of the problem can be obtained by contacting the corresponding
author.

##  III Result and discussion

In this section, we present normalized radiation pressure $F/F^{0}$, where
$F^{0}=\varepsilon E^{2}_{0}$, acting on an infinite dielectric (chiral) slab
by taking two linearly (circularly) polarized plane waves of same polarization
and vacuum wavelength $\lambda_{0}=1064{\rm nm}$. In all examples, we consider
refractive index of the slab $n_{s}=\sqrt{\varepsilon_{s}\mu_{s}}$, where
$\mu_{s}=1$ is the permeability of the slab, $d$ is the thickness of the slab,
and medium around the slab is water with relative refractive index
$n_{m}=1.332$.

Figure 4: (a) Normalized radiation pressure acting on a dielectric slab by two
plane waves impinging at $x=0$, with different incident angles:
$\theta=0\degree$ (dashed), $\theta=20\degree$ (red), $\theta=70\degree$
(black) as a function of $d/\lambda$. Normalized radiation pressure (b) as
function of $\theta$ and traversal position for fixed $d=1.5\lambda$, (c) as
function of $\theta$ and $d/\lambda$ for $x=0$, where, the color areas show
negative and white area shows positive radiation pressure.

First we consider a dielectric slab to elucidate the role of oblique
incidence, i.e., angle between the two impinging linearly polarized plane
waves, on the radiation pressure by using the analytically derived expression
given in Eq. (6). For the sake of consistency, one can get optical force
expression on a dielectric slab due to the normal incidence by setting
$\theta=0$. In this scenario, no interference occurs 1(b), and the incident
fields exert a positive acceleration on the center of mass of the slab by
virtue of the law of conservation of momentum. Thus, it manifests that the
radiation pressure by two incident normal plane waves is pushing for ${\rm
Im}(n_{s})\geq 0$ Mizrahi2010 ; Wang2015 ; Li2016 . Moreover, for better
illustration, we calculate the normalized radiation pressure as a function of
$d/\lambda$ for fixed real refractive index and the results are displayed in
Fig. 3: (a) ${\rm Re}{(n_{s})}=3.54$; and (b) ${\rm Re}{(n_{s})}=1.5$.
Whereas, in both the cases different values of imaginary part of the
refractive index have been used as: ${\rm Im}(n_{s})=0.05$ (red line), ${\rm
Im}(n_{s})=0$ (dashed line), and ${\rm Im}(n_{s})=-0.05$ (black line). It is
clear that for passive slab, radiation pressure is positive for ${\rm
Im}(n_{s})\geq 0$. On the other hand, pressure is negative when slab is active
Ali2021gain ; Ali2022acs ; Ali2021FIO  i. e. ${\rm Im}(n_{s})<0$ as
demonstrated in Fig. 3 (black lines). Since there is no interference and the
pressure is constant over the slab 1(b). Thus, the center of mass of the slab
undergoes to a net pushing force when slab is passive and pulling force when
slab is active. Another interesting feature is seen when we increase the
thickness of the slab, that is, the force is periodic for a lossless slab
because no photons are captured by the lossless slab, however, there is
periodically increasing positive pressure for the absorptive case, due to
absorptive losses. For active slab, the negative radiation pressure increases
by increasing the $d$ due to the fact that a thicker active slab emits more
light than it receives which allows for higher amplification of the forward
momentum, consequently increasing the negative radiation pressure.

It is worth mentioning that the slab with larger ${\rm Re}(n_{s})$ presents
larger reflection, therefore, it needs a thick active slab (i.e., large
emission requires to overcome the radiation losses) to revert the direction of
radiation pressure. For instance, in the Fig 3(a) ($n_{s}=3.54-0.05i$) the
radiation pressure changes sign at $d/\lambda=0.2$ and in Fig. 3(b)
($n_{s}=1.5-0.05i$) and the pressure becomes negative for small thickness.

The situation changes dramatically, when the two plane waves make an oblique
incidence at an angle $\theta>0\degree$, which lead to an interference effect
between the sources as shown in 1(c). As a result, the transmitted and
reflected lights give local maxima and minima along the x-direction and hence
the corresponding radiation pressure over the surface will provide periodic
locally positive and negative interference pattern 1(b). In order to observe
local positive and negative pressure on the slab we separate the both sources
by an angle $2\theta$ as shown in Fig. 1.

In Fig. 4, we consider two linearly polarized plane waves impinging on a
dielectric infinite slab of refractive index $n_{s}=1.5$. Fig. 4(a) presents
the radiation pressure as a function of $d/\lambda$ for different angles:
$\theta=0\degree$ (dashed), $\theta=20\degree$ (red), and $\theta=70\degree$
(black) for a fixed transversal position $x=0$. It is clear that for
$\theta=0$ and $\theta=20\degree$ radiation pressure is always positive which
means that scattered momentum is larger over the slab. However, for the large
angle the radiation pressure becomes negative which appears due to the
increase of forward momentum than the backward momentum and hence pressure is
locally negative. This locally positive and negative radiation pressure
pattern appears along the $x$-axis of infinite slab, which is shown in Fig.
4(b). It is clearly shown that for $\theta>30\degree$ the dielectric slab
undergoes to a local pulling and pushing force densities. The crossover
between the positive and negative radiation pressure is shown in Fig. 4(c),
where optical radiation pressure is calculated as a function of angle and
$d/\lambda$ for fixed $x=0$. It indicates that the slab at the point of
incidence experiences a positive radiation pressure $(F/F^{0}>0)$ for angle
$0\degree\leq\theta<30\degree$ and local negative radiation pressure
$(F/F^{0}<0)$ for $30\degree<\theta<90\degree$.

Figure 5: Normalized radiation pressure by two left circularly polarized plane
waves incident at $\theta=50\degree$ (green and red lines) acting on chiral
slab as a function of transversal position for different chirality parameter:
such as $\kappa=+0.3$ (red), $\kappa=-0.3$ (green). The dashed (emoji) line is
calculated for normal incident (i. e. $\theta=0\degree$) for chirality
parameter $\kappa=-0.3$ ($\kappa=0.3$). The thickness of the slab is fixed
$d=1.5\lambda$

.

For many applications, natural chiral materials such as proteins, organic
compounds, nucleic acids Corradini2007 , and chiral metamaterials and
metasurfaces are widely used in biological applications. In this context, if
the slab is made of such a chiral material, then a TE plane wave exerts an
identical radiation pressure on both chiral enantiomers with opposite
chirality parameter $\kappa$, as expected from symmetry Ali2021 . It is due to
fact that the TE plane wave reads the same refractive index regardless of
$\kappa$ and hence cannot create the specific positive and negative radiation
pressure intervals depending slightly on $\kappa$. To circumvent this issue we
come up with a new strategy, i.e., by shining two circularly polarized plane
waves on infinite chiral slab as shown in Fig 2. Now each polarization reads
distinct refractive index $n_{s}=\sqrt{\epsilon_{s}}+\sigma\kappa$, thus
exerted pressure on the chiral slab is obliged to be distinct.

Figure 6: Radiation pressure by two LCP plane waves impinging at an angle
$\theta=50\degree$ acting on dielectric slab as a function of $d/\lambda$ (a)
for chirality parameter $\kappa=-0.2$ (red) and $\kappa=0.2$ (green) (b)
$\kappa=-0.5$ (red) and $\kappa=0.5$ (green).

In Fig. 5 we demonstrate the normalized radiation pressure acting on a chiral
slab of thickness $d$ refractive index $n_{s}=1.5$ as a function of the
transversal axis for chirality parameter: $\kappa=-0.3$ and $\kappa=0.3$ at
different angles. First, we select normal incidence $\theta=0\degree$ and the
radiation pressure is calculated for $\kappa=0.3$ (dashed line) and
$\kappa=-0.3$ (emoji). It suggests that the radiation pressure is constant
along the tranversal direction, as expected (due to the vanishing of
interference) but it also manifests that the radiation pressure does not
depend on chirality for normal incidence, see dashed and emoji lines in Fig.
5. Thus, the electromagnetic field impinges at angle i. e. $\theta=50\degree$,
thereby, the interference comes into play. The chiral slab undergoes a local
positive and negative pressure, which explicitly depends on the chirality
parameters. For instance, left circularly polarized light exerts negative
(positive) radiation pressure on the chiral slab with
$\kappa=0.3(\kappa=-0.3)$ at origin (point of incidence), which propagates
along the transverse direction. The direction and magnitude of the radiation
pressure can be changed by changing the angle of the impinging plane waves and
the thickness of the slab.

In Fig. 6, we unveil the role of chirality by calculating the radiation
pressure as a function of $d/\lambda$ for different chirality parameters: Fig.
6(a) $\kappa=0.2$ (green) $\kappa=-0.2$ (red) Fig. 6(b) $\kappa=0.5$ (green)
$\kappa=-0.5$ (red) for fixed left circularly polarized light sources
impinging at $x=0$. For this particular configuration, chiral slab with
$\kappa>0$ (as the same handedness as of the light sources) undergoes to a
negative local radiation pressure when the handedness of chirality is reverted
(i.e. $\kappa<0$), then the pressure is positive. It can be associated with
the well know fact that a chiral medium favours the transmission of light if
the handed of the chiral is same as the handed of the polarization and on
contrast optimizes the reflection for opposite handedness.

This effect is elaborated in Fig 6(b) by taking large chirality parameters,
where radiation pressure corresponding to $\kappa=-0.5$ presents a maximum
positive pressure because the maximum of the incident signal gets reflected
due to opposite handedness and the radiation pressure corresponding to
$\kappa=0.5$ (green line) provides a maximum negative pressure, as expected by
the symmetry consideration. This pressure can be reversed by changing the
polarization of the incident beam. Indeed, the radiation pressure is invariant
by changing both $\kappa$ and $\sigma$ in the same way.

##  IV Conclusions

In conclusion, we have developed, for the first time, a theoretical framework
to calculate the local radiation pressure on an infinite dielectric (chiral)
slab illuminated by a superposition of linearly (circularly) polarized plane
waves. In this prospective, we have derived explicit analytical results for
the optical radiation pressure acting on an infinite slab by following two
approaches. In the first approach, we have solved the standard Maxwell stress
tensor to calculate the total optical radiation pressure by incorporating the
total fields on each side of the slab. Consequently, the slab experiences a
total local positive or negative radiation pressure. In the second approach,
we have calculated the Maxwell stress tensor at each interface of the slab
separately and then obtained the explicit expressions for the radiation
pressure on respective interfaces. However, subsequently, the superposition of
the radiation pressures of the both interfaces provides the total radiation
pressure on the slab. Moreover, we have not only verified that the final
results obtained by the both approaches are in good agreement but we also have
shown, by considering the limiting cases (for instance, by taking the incident
angle equal to zero), that our analytical results perfectly agree with the
well-known results of the radiation pressure acting on a slab due to normal
impinging field. In addition to the analytical work, we have considered
realistic material parameters to numerically validate our findings. These
results have shown that the local radiation pressure strongly depends on the
incident angle, polarization of the light source and the nature of the
material. For instance, for a normal incidence of the field, the radiation
pressure is constant all over the slab, whereas, for an oblique incidence it
varies locally along the transversal direction, which occurs due to the
interference effects. In the case of chiral slab, chirality parameter also
plays a key role to govern the local minima and maxima of the radiation
pressure, which are indeed consistent with the momentum and energy
conservation laws. At last but not the least, we hope that this approach paves
the path for calculating the radiation pressure that can be used, in
conjunction with numerical simulations, to yield the distribution of fields,
as well as the resulting interference patterns and corresponding local
extension and relaxation of soft matter, cell membranes in diverse biological
systems of practical interest.

## Appendix

##  Appendix A Force at the interfaces

In this appendix we present the principal steps to obtain the radiation
pressure due to two non-collinear plane waves on each interfaces of the slab
as shown in Fig. 1.

In order to calculate the radiation pressure $\mathbf{F}(z=0,x,\theta)$ at the
interface $z=0$, we take the flux through the Gaussian surface at $z=0$ and
the resultant $\mathbf{F}(z=0,x,\theta)$ is expression in terms of the Maxwell
stress tensor as

$\mathbf{F}(z=0,x,\theta)=\frac{1}{2}\Re\\{\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=0^{-})-\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=0^{+})\\}.$
(18)

In this case the electromagnetic fields inside the slab can be expressed by
two types of plane waves propagating with positive and negative $z$ components
of wave vectors

$\mathbf{E}_{w,k_{z}>0}=2E_{w}^{(1)}cos(k_{s}x\sin\theta_{s})e^{ik_{s}z\cos\theta_{s}}e^{-i\omega
t}\mathbf{\hat{y}}$ (19)

and

$\mathbf{E}_{w,k_{z}<0}=2E_{w}^{(2)}cos(k_{s}x\sin\theta_{s})e^{-ik_{s}z\cos\theta_{s}}e^{-i\omega
t}\mathbf{\hat{y}},$ (20)

respectively, where $k_{s}$ is the wave vector inside slab. In addition,

$E_{w}^{(1)}=\frac{E_{0}t_{m,s}}{1-r_{s,m}^{2}\,e^{i2\beta}}$ (21)

and

$E_{w}^{(2)}=\frac{E_{0}t_{m,s}r_{s,m}e^{i\beta}}{1-r_{s,m}^{2}\,e^{i2\beta}}.$
(22)

The total electric field is given by
$\mathbf{E}_{total,w}=\mathbf{E}_{w,k_{z}>0}+\mathbf{E}_{w,k_{z}<0}$.
Furthermore, the total electric fields at the both sides of the Gaussian
interface $z=0$ are evaluated,
$\mathbf{E}_{total}(z=0^{-})=\mathbf{E}_{in,total}(z=0^{-})+\mathbf{E}_{r,total}(z=0^{-})$
and $\mathbf{E}_{total}(z=0^{+})=\mathbf{E}_{total,w}(z=0^{+})$. By setting
these fields along with Eqs. 1-2 and 19-20 in Eq.5 we obtain the radiation
pressure acting on at the interface $z=0$:

$\displaystyle\mathbf{F}(z=0,x,\theta)=\varepsilon
E_{0}^{2}\biggl{[}\biggl{(}|1+r_{TE}|^{2}+|1-r_{TE}|^{2}\cos^{2}\theta\biggr{)}$
$\displaystyle\cos^{2}(kx\sin\theta)-|1+r_{TE}|^{2}\,\sin^{2}\theta\,\sin^{2}(kx\sin\theta)\biggr{]}\mathbf{\hat{z}}$
$\displaystyle-\,\varepsilon_{s}E_{0}^{2}\biggl{[}\biggl{(}|t_{w}+r_{w}e^{i\beta}|^{2}+|t_{w}-r_{w}e^{i\beta}|^{2}\cos^{2}\theta_{s}\biggr{)}$
$\displaystyle\cos^{2}(k_{s}x\sin\theta_{s})-|t_{w}+r_{w}e^{i\beta}|^{2}\,\sin^{2}\theta_{s}\,\sin^{2}(k_{s}x\sin\theta_{s})\biggr{]}\mathbf{\hat{z}},$
(23)

where

$t_{w}=\frac{t_{m,s}}{1-r_{s,m}^{2}e^{i2\beta}}$ (24)

and

$r_{w}=\frac{t_{m,s}r_{s,m}}{1-r_{s,m}^{2}e^{i2\beta}}.$ (25)

Similarly, the radiation pressure $\mathbf{F}(z=d,x,\theta)$ at the interface
$z=d$ is given by

$\mathbf{F}(z=d,x,\theta)=\frac{1}{2}\Re\\{\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=d^{-})-\mathbf{\hat{z}}\cdot\overleftrightarrow{T}(z=d^{+})\\}.$
(26)

The total electric fields in the both sides of the gaussian surface are given
by $\mathbf{E}_{total}(z=d^{-})=\mathbf{E}_{total,w}(z=d^{-})$ and
$\mathbf{E}_{total}(z=d^{+})=\mathbf{E}_{t,total}(z=d^{+})$, and the radiation
pressure is obtained setting the Eqs. 3 and 19-20 in Eq.5

$\displaystyle\mathbf{F}(z=d,x,\theta)=-\varepsilon
E_{0}^{2}|r_{TE}|^{2}\biggl{[}\biggl{(}1+\cos^{2}\theta\biggr{)}\,\cos^{2}(kx\sin\theta)$
$\displaystyle-\sin^{2}\theta\,\sin^{2}(kx\sin\theta)\biggr{]}\mathbf{\hat{z}}+\,\varepsilon_{s}E_{0}^{2}\biggl{[}\biggl{(}|t_{w}e^{i\beta}+r_{w}|^{2}$
$\displaystyle+|t_{w}e^{i\beta}-r_{w}|^{2}\cos^{2}\theta_{s}\biggr{)}\,\cos^{2}(k_{s}x\sin\theta_{s})-|t_{w}e^{i\beta}+r_{w}|^{2}$
$\displaystyle\sin^{2}\theta_{s}\sin^{2}(k_{s}x\sin\theta_{s})\biggr{]}\mathbf{\hat{z}}.$
(27)

Finally, the total radiation pressure acting the slab is written as

$\mathbf{F}=\mathbf{F}(x,\theta)=\mathbf{F}(z=0,x,\theta)+\mathbf{F}(z=d,x,\theta)$
(28)

It is clear from Eqs. (23) and (27) that the net momentum flux inside the slab
is conserved, and the total radiation pressure on the surface is given in Eq.
(28) which is the same as we have presented in Eq. (6). Thus, by analysing the
Eq. (28), it clear that one can not only recover Eq. (6) for the case of
losses slab but also recover Eq. (12), for normal incidence even considering
the case when imaginary part of the refractive index is different than zero
$Im(n_{s})\neq 0$.

## Funding

We thank G. Wiederhecker, F. A. Pinheiro, F. S. S. Rosa and P. A. M. Neto for
inspiring discussion. This work is partially supported by Fundacão de Amparo a
Pesquisa do Estado de São Paulo (FAPESP) (2020/03131-2). The authors declare
no conflicts of interest.

## References

  * (1) A. Ashkin, Acceleration and trapping of particles by radiation pressure, Phys. Rev. Lett. 24, 156 (1970).
  * (2) A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Observation of a single-beam gradient force optical trap for dielectric particles, Opt. Lett. 11, 288 (1986).
  * (3) J. Gieseler, J. R. Gomez-Solano, et al., Optical Tweezers: A Comprehensive Tutorial from Calibration to Applications, Adv. Opt. Photon. 13, 74-241 (2021).
  * (4) P. Polimeno, A. Magazzú, M. A. Iatí, et al.,Optical tweezers and their applications, J. Quant. Spectrosc. Radiat. Transf. 218, 131-150 (2018).
  * (5) A. Ashkin, Optical trapping and manipulation of neutral particles using lasers: A reprint volume with commentaries (World Scientific, Singapore, 2006).
  * (6) M. E. J. Friese, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, Optical alignment and spinning of laser-trapped microscopic particles, Nature 394, 348 (1998).
  * (7) D. G. Grier, A revolution in optical manipulation, Nature 424, 810 (2003).
  * (8) J. Lin and Y. Q. Li, Optical trapping and rotation of airborne absorbing particles with a single focused laser beam, Appl. Phys. Lett. 104, 101909 (2014).
  * (9) R. Ali, F. A. Pinheiro, R. S. Dutra, and P. A. Maia Neto, Tailoring optical pulling forces with composite microspheres, Phys. Rev. A 102, 023514 (2020).
  * (10) R. Ali, R. S. Dutra, F. A. Pinheiro, and P. A. Maia Neto, Enantioselection and chiral sorting of single microspheres using optical pulling forces , Opt. Lett. 46, 1640 (2021).
  * (11) J. Chen, J. Ng, Z. F. Lin, and C. T. Chan, Optical pulling force, Nat. Photonics 5, 531 (2011).
  * (12) X. Li, J. Chen, Z. Lin, J. Ng. Optical pulling at macroscopic distances, Sci. Adv. 5, 7814 (2019).
  * (13) V. Shvedov, A. R. Davoyan, N. Engheta and W. Krolikowski, A long-range polarization-controlled optical tractor beam, Nat. Photonics  8, 846 (2014).
  * (14) S. Sukhov and A. Dogariu, Negative nonconservative forces: optical tractor beams for arbitrary objects, Phys. Rev. Lett. 107, 203602 (2011).
  * (15) R. Ali, F. A. Pinheiro, F. S. S. Rosa, R. S. Dutra, P. A. Maia Neto, Optimizing optical tweezing with directional scattering in composite microspheres, Phys. Rev. A 98 , 053804 (2018).
  * (16) S. H. Simpson, and and S. Hanna, Optical trapping of spheroidal particles in Gaussian beams , J. Opt. Soc. Am. A 24, 430-443 (2007).
  * (17) B. J. Wiley, Y. Chen, J. M. McLellan, Y. Xiong, Z.-Y. Li, D. Ginger, and Y. Xia, Synthesis and Optical Properties of Silver Nanobars and Nanorice,  Nano Letters, 7, 1032-1036 (2007).
  * (18) R. C. Gauthier, M. Ashman, and C. P. Grover, Experimental confirmation of the optical-trapping properties of cylindrical objects, Appl. Opt. 38, 4861-4869 (1999).
  * (19) R. Ali, F. A. Pinheiro, F. S. S. Rosa, R. S. Dutra, P. A. Maia Neto,  Enantioselective manipulation of chiral nanoparticles using optical tweezers, Nanoscale 12, 5031 (2020).
  * (20) R. Ali, F. A. Pinheiro, R. S. Dutra, F. S. S. Rosa, and P. A. Maia Neto, Probing the optical chiral response of single nanoparticles with optical tweezers, J. Opt. Soc. Am. B 37, 2796-2803 (2020).
  * (21) A. Canaguier-Durand, J. A. Hutchison, C. Genet, and T. W. Ebbesen, Mechanical separation of chiral dipoles by chiral light, New J. Phys. 15, 123037 (2013).
  * (22) K. Ding, J. Ng, L. Zhou, and C. T. Chan,  Realization of optical pulling forces using chirality, Phys. Rev. A 89, 063825 (2014).
  * (23) M. Mansuripur, Radiation pressure and the linear momentum of the electromagnetic field, Opt. Express 12, 5375-5401 (2004).
  * (24) R. Loudon and S. M. Barnett, Theory of the radiation pressure on dielectric slabs, prisms and single surfaces, Opt. Express 14, 11855-11869 (2006).
  * (25) R. Alaee, J. Christensen, and M. Kadic, Optical Pulling and Pushing Forces in Bilayer PT -Symmetric Structures, Phys. Rev. Appl. 9, 014007 (2018)
  * (26) , R. Corradini, S. Sforza, T. Tedeschi, and R. Marchelli, Chirality as a tool in nucleic acid recognition: Principles and relevance in biotechnology and in medicinal chemistry, Chirality 19, 269, (2007)
  * (27) A. Basiri, X. Chen, J. Bai, et al. Nature-inspired chiral metasurfaces for circular polarization detection and full-Stokes polarimetric measurements, Light Sci. Appl. 8, 78 (2019).
  * (28) S. Lepeshov, and A. Krasnok, Virtual optical pulling force, Optica, 7, 1024-1030 (2020).
  * (29) A. Mizrahi, and Y. Fainman, Negative radiation pressure on gain medium structures, Opt. Lett. 35, 3405-3407 (2010).
  * (30) M. Wang, H. Li, D. Gao, L. Gao, J. Xu, and Cheng-Wei Qiu, Radiation pressure of active dispersive chiral slabs, Opt. Express 23, 16546-16553 (2015).
  * (31) G. Li, M. Wang, H. Li, M. Yu, Y. Dong, and Jun Xu, Wave propagation and Lorentz force density in gain chiral structures, Opt. Mater. Express 6, 388-395 (2016).
  * (32) R. Ali, R. S. Dutra, F. A. Pinheiro and P. A. Maia Neto, Gain-assisted optical tweezing of plasmonic and large refractive index microspheres, J. Opt. 23, 115004, (2021).
  * (33) M. Born, and E. Wolf. Principles of optics: electromagnetic theory of propagation, interference and diffraction of light, Elsevier, (2013).
  * (34) A. B. Kemp, M. T. Grzegorczyk, and J. A. Kong, Ab initio study of the radiation pressure on dielectric and magnetic media, Opt. Express 13, 9280-9291 (2005).
  * (35) R. N. C. Pfeifer, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, Momentum of an electromagnetic wave in dielectric media, Rev. Mod. Phys. 79, 1197 (2007).
  * (36) S. B. Wang and C.T. Chan, Lateral optical force on chiral particles near a surface, Nat. Commun. 5, 3307 (2014).
  * (37) C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, (Wiley, New York, 1983).
  * (38) R. Ali, R. S. Dutra, P. A. Pinheiro, F. S. S. Rosa, and P. M. A. Neto, Theory of optical tweezing of dielectric microspheres in chiral host media and its applications, Sci. Rep. 10, 16481 (2020).
  * (39) R. Ali, T. P. M. Alegre, F. A. Pinheiro, and G. S. Wiederhecker, ”Enantioselective optical forces of gain functionalized core-shell chiral nanoparticles,” in Frontiers in Optics + Laser Science 2021, C. Mazzali, T. (T.-C.) Poon, R. Averitt, and R. Kaindl, eds., Technical Digest Series (Optica Publishing Group, 2021), paper FTh4E.2.
  * (40) R. Ali, Tunable anomalous scattering and negative asymmetry parameter in a gain-functionalized low refractive index sphere, ACS Omega 7, 2170–2176 (2022).
  * (41) S. Bassiri, C. H. Papas, and N. Engheta,  Electromagnetic wave propagation through a dielectric chiral interface and through a chiral slab, J. Opt. Soc. Am. A,  5, 1450-1459, (1988).


# Limit cycles of linear vector fields on
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$

###### Abstract

It is well known that linear vector fields defined in $\mathbb{R}^{n}$ can not
have limit cycles, but this is not the case for linear vector fields defined
in other manifolds. We study the existence of limit cycles bifurcating from a
continuum of periodic orbits of linear vector fields on manifolds of the form
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$ when such vector fields are
perturbed inside the class of all linear vector fields. The study is done
using the averaging theory. We also present an open problem concerning the
maximum number of limit cycles of linear vector fields on
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$.

Clara Cufí-Cabré and Jaume Llibre

Departament de Matemàtiques

Universitat Autònoma de Barcelona (UAB)

08193 Bellaterra, Barcelona, Spain

<EMAIL_ADDRESS><EMAIL_ADDRESS>

_2020 Mathematics Subject Classification:_ 34C29, 34C25, 34A30.
_Key words:_ limit cycle, periodic orbit, isochronous center, averaging
method.

##  1 Introduction and statement of the main results

The study of periodic orbits of differential systems play an important role in
the qualitative theory of ordinary differential equations and their
applications. A _limit cycle_ is defined as a periodic orbit of a differential
system which is isolated in the set of all periodic orbits of the system.
There are many works concerning the study of limit cycles and their
applications (see for instance [3, 5, 6, 8] and the references quoted
therein).

It is well know that linear vector fields in $\mathbb{R}^{n}$ can not have
limit cycles, but this is not the case if one considers linear vector fields
in other manifolds different from $\mathbb{R}^{n}$. The objective of this
paper is to study the existence of limit cycles of linear vector fields
defined on the manifolds $(\mathbb{S}^{2})^{n}\times\mathbb{R}^{n}$.

The problem of studying limit cycles of linear vector fields on manifolds
different form $\mathbb{R}^{n}$ was already treated in [14], where the authors
consider linear vector fields on $\mathbb{S}^{m}\times\mathbb{R}^{n}$, and
they conjecture that such vector fields may have at most one limit cycle.

Linear autonomous differential systems, namely, systems of the form
$\dot{x}=Ax+b$, where $A$ is a $n\times n$ real matrix and $b$ is a vector in
$\mathbb{R}^{n}$, are the easiest systems to study because their solutions can
be completely determined (see [1, 18]), but still they play an important role
in the theory of differential systems. Thus when a nonlinear differential
system has a hyperbolic equilibrium point, the dynamics around that point is
determined by the linearization of the vector field at that point (Hartman-
Grossman theorem, see [7]).

Also linear vector fields having invariant subspaces of periodic orbits can be
perturbed inside a concrete class of nonlinear differential systems to obtain
limit cycles of these nonlinear systems bifurcating from the periodic orbits
of the linear system (see [4, 11, 12, 13]).

Moreover, linear differential systems of the form $\dot{x}=Ax+Bu$, where $x$
are the state variables and $u$ is the control input, are applied in control
theory for the modeling of hybrid systems (see [9, 10]).

The preceeding exemples show some of the importance of linear differential
systems. In this paper we show that linear differential systems can have limit
cycles when the manifold where they are defined is different from
$\mathbb{R}^{n}$, and we consider the question of how many limit cycles can
have at most a linear vector field depending on the manifold where it is
defined.

Let $M$ be a smooth connected manifold of dimension $n$, and let $TM$ be its
tangent bundle. A _vector field_ on $M$ is a map $X:M\to TM$ such that
$X(x)\in T_{x}M$, where $T_{x}M$ is the tangent space of $M$ at the point $x$.

A _linear vector field_ in $\mathbb{R}^{n}$ is a vector field of the form
$X(x)=Ax+b$, with $x,b\in\mathbb{R}^{n}$ and where $A$ is a $n\times n$ real
matrix. As it is well known linear vector fields on $\mathbb{R}^{n}$ either do
not have periodic orbits or their periodic orbits form a continuum, and
therefore they do not have limit cycles.

In this paper we consider linear vector fields on some manifolds of the form
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$, where $\mathbb{S}^{2}$ denotes the
unit two-dimensional sphere. Here the sphere $\mathbb{S}^{2}$ is parameterized
by the coordinates $(\theta,\varphi)$, where $\theta\in[-\pi,\pi)$ denotes the
azimuth angle and $\varphi\in[-\pi/2,\pi/2]$ is the polar angle. Hence the
curve $\\{\varphi=0\\}$ is the equator of the sphere.

Let $(\theta_{1},\varphi_{1},\dots\theta_{m},\varphi_{m},x_{1},\cdots,x_{n})$
denote the coordinates of the space
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$. Then we say that a vector field
$X$ is _linear_ on $M=(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$ if the
expression of $X$ in the coordinates
$z=(\theta_{1},\varphi_{1},\dots\theta_{m},\varphi_{m},x_{1},\cdots,x_{n})\in
M$ is of the form $X(z)=Az+b$, with $b\in M$ and where $A$ is a
$(2m+n)\times(2m+n)$ real matrix.

A simple example in which a linear differential system on the manifold
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$ has a limit cycle is the following.
Take $m=1,~\,n=0$ and consider the linear system on the sphere
$\mathbb{S}^{2}$ given by

$\dot{\theta}=1,\qquad\dot{\varphi}=\varphi,$

for $\theta\in[-\pi,\pi)$ and $\varphi\in(-\pi/2,\pi/2)$, and

$\dot{\theta}=0,\qquad\dot{\varphi}=0$

for $\varphi=\pm\pi/2$. Then, clearly the equator of the sphere
$\\{\varphi=0\\}$ is the only periodic orbit of the system, and therefore it
is a limit cycle.

In this paper we consider generic linear perturbations of some linear vector
fields on three different manifolds of the form
$(\mathbb{S}^{2})^{m}\times\mathbb{R}^{n}$, and we study whether those
families of linear differential systems can have limit cycles.

Let $M=\mathbb{S}^{2}\times\mathbb{R}$ and consider the linear differential
system in $M$ given by

$\dot{\theta}=1,\quad\dot{\varphi}=0,\quad\dot{r}=r-1,$ (1.1)

for $r\in\mathbb{R}$, $\theta\in[-\pi,\pi)$ and $\varphi\in(-\pi/2,\pi/2)$,
and with $\dot{\theta}=0$ on the straight lines $R_{1}=\\{\varphi=-\pi/2\\}$
and $R_{2}=\\{\varphi=\pi/2\\}$.

The solution of system (1.1) is given by

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad
r(t)=(r_{0}-1)e^{t}+1.$

Thus the sphere $\\{r=1\\}$ is an invariant manifold with two equilibrium
points at the north and the south poles, and is foliated by periodic orbits of
period $2\pi$, corresponding to the parallels of the sphere, except at the
poles. Moreover the straight lines $R_{1}$ and $R_{2}$ are invariant.

First we shall study the bifurcation of limit cycles when we perturb system
(1.1) inside the class of all linear differential systems, and we shall see
that one of the periodic orbits contained in the sphere $\\{r=1\\}$ may
bifurcate to a limit cycle under certain hypotheses.

We consider the class of differential systems

$\displaystyle\begin{split}&\dot{\theta}=1+\varepsilon(a_{0}+a_{1}\theta+a_{2}\varphi+a_{3}r),\\\
&\dot{\varphi}=\varepsilon(b_{0}+b_{1}\theta+b_{2}\varphi+b_{3}r),\\\
&\dot{r}=r-1+\varepsilon(c_{0}+c_{1}\theta+c_{2}\varphi+c_{3}r).\end{split}$
(1.2)

where $a_{i}$, $b_{i}$ and $c_{i}$, for $i=0,\dots,3$ are real numbers and
with $\varepsilon>0$ being a small parameter. Note that this is the more
general linear perturbation of system (1.1). For the class of systems
(LABEL:pertorb-lineal-2) we have the following result.

######  Theorem 1.

For sufficiently small $\varepsilon>0$ the linear differential system
(LABEL:pertorb-lineal-2) has a limit cycle bifurcating from a periodic orbit
of system (1.1) provided that $a_{1}b_{2}-a_{2}b_{1}\neq 0$. Moreover this
limit cycle bifurcates from the periodic orbit of system (1.1) parameterized
by $(\theta(t),\varphi(t),r(t))=(\theta_{0}+t,\varphi_{0},1)$, with

$\displaystyle\theta_{0}$
$\displaystyle=\frac{a_{2}(b_{0}+b_{3}+b_{1}\pi)-b_{2}(a_{0}+a_{3}+a_{1}\pi)}{a_{1}b_{2}-a_{2}b_{1}},$
$\displaystyle\varphi_{0}$
$\displaystyle=\frac{b_{1}(a_{0}+a_{3}+a_{1}\pi)-a_{1}(b_{0}+b_{3}+b_{2}\pi)}{a_{1}b_{2}-a_{2}b_{1}}.$

Theorem 1 is proved in section 3.

We remark that the existence of the limit cycle for system (LABEL:pertorb-
lineal-2) does not depend on the perturbation of the $\dot{r}$ equation.

As an example of the previous result, consider the system

$\dot{\theta}=1+\varepsilon a\varphi,\quad\dot{\varphi}=\varepsilon
b\theta,\quad\dot{r}=r-1,$ (1.3)

with $a,b\in\mathbb{R}$ and $\varepsilon>0$. In this case the sphere
$\\{r=1\\}$ is still an invariant manifold. Appliying Theorem 1 with
$a_{2}=a,\,b_{1}=b$ and the rest of the coefficients of the perturbation being
zero, we find that system (1.3) has a limit cycle bifurcating form the
periodic orbit of system (1.1) parameterized by
$(\theta(t),\varphi(t),r(t))=(-\pi+t,0,1)$. That is, there is a limit cycle
bifurcating from the periodic orbit corresponding to the equator of the sphere
$\\{r=1\\}$ of system (1.1). Moreover this limit cycle is still contained in
the sphere $\\{r=1\\}$.

Next we consider linear differential systems defined on higher dimensional
manifolds. We take $M=\mathbb{S}^{2}\times\mathbb{S}^{2}\times\mathbb{R}$ and

$\dot{\theta}=1,\quad\dot{\varphi}=0,\quad\dot{\nu}=1,\quad\dot{\phi}=0,\quad\dot{r}=r-1,$
(1.4)

for $(\theta,\varphi,\nu,\phi,r)\in M$, with $\theta,\nu\in[-\pi,\pi)$ and
$\varphi,\phi\in(-\pi/2,\pi/2)$, and with $\dot{\theta}=0$ when
$\varphi=\pm\pi/2$ and $\dot{\nu}=0$ when $\phi=\pm\pi/2.$

The general solution of system (1.4) is

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad\nu(t)=\nu_{0}+t,\quad\phi(t)=\phi_{0},\quad
r(t)=(r_{0}-1)e^{t}+1,$

and thus the product of spheres $\\{r=1\\}\cong(\mathbb{S}^{2})^{2}$ is an
invariant manifold foliated by periodic orbits of period $2\pi$, except for
the four points $\\{r=1,\varphi=\pm\pi/2,\phi=\pm\pi/2\\}$, which are
equilibrium points.

We consider the most general perturbation of the differential system (1.4)
inside the class of all linear differential systems, namely

$\displaystyle\begin{split}&\dot{\theta}=1+\varepsilon(a_{0}+a_{1}\theta+a_{2}\varphi+a_{3}\nu+a_{4}\phi+a_{5}r),\\\
&\dot{\varphi}=\varepsilon(b_{0}+b_{1}\theta+b_{2}\varphi+b_{3}\nu+b_{4}\phi+b_{5}r)\\\
&\dot{\nu}=1+\varepsilon(c_{0}+c_{1}\theta+c_{2}\varphi+c_{3}\nu+c_{4}\phi+c_{5}r),\\\
&\dot{\phi}=\varepsilon(d_{0}+d_{1}\theta+d_{2}\varphi+d_{3}\nu+d_{4}\phi+d_{5}r)\\\
&\dot{r}=r-1+\varepsilon(e_{0}+e_{1}\theta+e_{2}\varphi+e_{3}\nu+e_{4}\phi+e_{5}r),\end{split}$
(1.5)

with $a_{i},b_{i},c_{i},d_{i},e_{i}\in\mathbb{R}$ for $i=0,\ldots,5$, and with
$\varepsilon>0$ being a small parameter. In the following result we give
sufficient conditions on the coefficients of system (LABEL:pertorb-lineal-
high) in order that there is a limit cycle bifurcating from a periodic orbit
of the corresponding unperturbed system.

######  Theorem 2.

For sufficiently small $\varepsilon>0$ the differential system (LABEL:pertorb-
lineal-high) has a limit cycle bifurcating from a periodic orbit of system
(1.4) provided that

$\det\begin{pmatrix}a_{1}&a_{2}&a_{3}&a_{4}\\\ b_{1}&b_{2}&b_{3}&b_{4}\\\
c_{1}&c_{2}&c_{3}&c_{4}\\\ d_{1}&d_{2}&d_{3}&d_{4}\end{pmatrix}\neq 0.$

Moreover this limit cycle bifurcates from the periodic orbit of system (1.4)
parameterized by
$(\theta(t),\varphi(t),\nu(t),\phi(t),r(t))=(\theta_{0}+t,\varphi_{0},\nu_{0}+t,\phi_{0},1)$,
where $(\theta_{0},\,\varphi_{0},\,\nu_{0},\,\phi_{0})$ is the unique solution
of the linear system

$\displaystyle\begin{split}&a_{1}\theta_{0}+a_{2}\varphi_{0}+a_{3}\nu_{0}+a_{4}\phi_{0}=-a_{0}-a_{1}\pi-
a_{3}\pi-a_{5},\\\
&b_{1}\theta_{0}+b_{2}\varphi_{0}+b_{3}\nu_{0}+b_{4}\phi_{0}=-b_{0}-b_{1}\pi-
b_{3}\pi-b_{5},\\\
&c_{1}\theta_{0}+c_{2}\varphi_{0}+c_{3}\nu_{0}+c_{4}\phi_{0}=-c_{0}-c_{1}\pi-
c_{3}\pi-c_{5},\\\
&d_{1}\theta_{0}+d_{2}\varphi_{0}+d_{3}\nu_{0}+d_{4}\phi_{0}=-d_{0}-d_{1}\pi-
d_{3}\pi-d_{5}.\\\ \end{split}$

Theorem 2 is proved in section 4.

Finally we consider the linear differential system defined in
$M=\mathbb{R}^{2}\times\mathbb{S}^{2}$, for
$(x,y,\theta,\varphi)\in\mathbb{R}^{2}\times\mathbb{S}^{2}$, with
$\theta\in[-\pi,\pi)$ and $\varphi\in(-\pi/2,\pi/2)$, given by

$\dot{x}=-y,\quad\dot{y}=x,\quad\dot{\theta}=1,\quad\dot{\varphi}=0,$ (1.6)

and with $\dot{\theta}=0$ in the planes $P_{1}=\\{\varphi=-\pi/2\\}$ and
$P_{2}=\\{\varphi=\pi/2\\}$, which are invariant. The general solution of
system (1.6) is

$x(t)=x_{0}\cos t-y_{0}\sin t,\quad y(t)=x_{0}\sin t+y_{0}\cos
t,\quad\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},$

and therefore the whole phase space is filled by periodic orbits of period
$2\pi$, except for the two equilibrium points
$(x,y,\theta,\varphi)=(0,0,\theta,-\pi/2)$ and
$(x,y,\theta,\phi)=(0,0,\theta,\pi/2)$.

We consider the most general linear perturbation of system (1.6) and we study
the existence of limit cycles bifurcating from the periodic orbits of system
(1.6).

Let

$\displaystyle\begin{split}&\dot{x}=-y+\varepsilon(a_{0}+a_{1}x+a_{2}y+a_{3}\theta+a_{4}\varphi),\\\
&\dot{y}=x+\varepsilon(b_{0}+b_{1}x+b_{2}y+b_{3}\theta+b_{4}\varphi),\\\
&\dot{\theta}=1+\varepsilon(c_{0}+c_{1}x+c_{2}y+c_{3}\theta+c_{4}\varphi),\\\
&\dot{\varphi}=\varepsilon(d_{0}+d_{1}x+d_{2}y+d_{3}\theta+d_{4}\varphi),\\\
\end{split}$ (1.7)

be the perturbed system, with $a_{i},b_{i},c_{i},d_{i}\in\mathbb{R}$ for
$i=0,\ldots,4$, and where $\varepsilon>0$ is a small parameter. For this
linear differential system we have the following result.

######  Theorem 3.

For sufficiently small $\varepsilon>0$ the linear differential system
(LABEL:pertorb-lineal-r2s2) has a limit cycle bifurcating from a periodic
orbit of system (1.6) provided that

$\det\begin{pmatrix}b_{2}+a_{1}&a_{2}-b_{1}\\\
b_{1}-a_{2}&b_{2}+a_{1}\end{pmatrix}\neq
0\quad\text{and}\quad\det\begin{pmatrix}c_{3}&c_{4}\\\
d_{3}&d_{4}\end{pmatrix}\neq 0.$

Moreover this limit cycle bifurcates form the periodic orbit of system (1.6)
passing through the point $(x_{0},y_{0},\theta_{0},\varphi_{0})$ where

$\displaystyle x_{0}$
$\displaystyle=\frac{(2b_{2}+2a_{1})b_{3}-2a_{3}b_{1}+2a_{2}a_{3}}{b_{2}^{2}+b_{1}^{2}+a_{2}^{2}+a_{1}^{2}+2a_{1}b_{2}-2a_{2}b_{1}},$
$\displaystyle y_{0}$
$\displaystyle=-\frac{(2b_{1}-2a_{2})b_{3}+2a_{3}b_{2}+2a_{1}a_{3}}{b_{2}^{2}+b_{1}^{2}+a_{2}^{2}+a_{1}^{2}+2a_{1}b_{2}-2a_{2}b_{1}},$
$\displaystyle\theta_{0}$ $\displaystyle=-\frac{(\pi c_{3}+c_{0})d_{4}-\pi
c_{4}d_{3}-c_{4}d_{0}}{c_{3}d_{4}-c_{4}d_{3}},$ $\displaystyle\varphi_{0}$
$\displaystyle=\frac{c_{0}d_{3}-c_{3}d_{0}}{c_{3}d_{4}-c_{4}d_{3}}.$

Theorem 3 is proved in section 5.

As an example consider the system

$\dot{x}=-y+\varepsilon ay,\quad\dot{y}=x+\varepsilon
bx,\quad\dot{\theta}=1+\varepsilon c\varphi,\quad\dot{\varphi}=\varepsilon
d\theta,$ (1.8)

with $a,b,c,d\in\mathbb{R}$, and $\varepsilon>0$. Appliying Theorem 3 with
$a_{2}=a,\,b_{1}=b,\,c_{4}=c,\,d_{3}=d$ and the rest of the coefficients of
the perturbation being zero, we obtain that system (1.8) has a limit cycle
bifurcating form the periodic orbit of system (1.6) passing through the point
$(x_{0},y_{0},\theta_{0},\varphi_{0})=(0,0,-\pi,0)$, provided that
$(a-b)cd\neq 0$. That is, here the limit cycle bifurcates from the periodic
orbit corresponding to the equator of the invariant sphere $\\{x=y=0\\}$ of
system (1.6).

The key tool that we use for proving Theorems 1 \- 3 is the averaging theory.
For a general introduction to this theory, see the books [17, 19]. As one can
see in the proofs of Theorems 1 \- 3, our method based on the averaging theory
can produce at most one limit cycle for the studied systems. Therefore the
following open question is natural.

Open question. Let $m$ and $n$ be two non-negative integers. Is it true that a
linear vector field on the manifold $(\mathbb{S}^{m})^{m}\times\mathbb{R}^{n}$
can have at most one limit cycle?

A similar open question was stated in [14] concerning linear vector fields on
the manifold $(\mathbb{S}^{1})^{m}\times\mathbb{R}^{n}$.

##  2 Basic results on the averaging theory

In this section we state some basic results from the averaging theory that
will be used to prove the main results of the paper.

Let $M$ be a smooth connected manifold of dimension $n$, and let
$F_{0},\,F_{1}:\mathbb{R}\times M\to\mathbb{R}^{n}$ and
$F_{2}:\mathbb{R}\times M\times[0,\varepsilon_{0})\to\mathbb{R}^{n}$ be
$C^{2}$, $T$-periodic functions. Given the differential system

$\dot{x}(t)=F_{0}(t,x),$ (2.1)

we consider a perturbation of this system of the form

$\dot{x}(t)=F_{0}(t,x)+\varepsilon
F_{1}(t,x)+\varepsilon^{2}F_{2}(t,x,\varepsilon).$ (2.2)

The objective is to study the bifurcation of $T$-periodic solutions of system
(2.2) for $\varepsilon>0$ small enough. A solution to this problem is given by
the averaging theory.

We assume that there exists $k\leq n$ such that $M=M_{k}\times M_{n-k}$, where
$M_{k}$ is a manifold of dimension $k$ and $M_{n-k}$ is a manifold of
dimension $n-k$, and that the unperturbed system, namely system (2.1),
contains an open set, $V\subseteq M_{k}$, such that $\overline{V}$ is filled
with periodic solutions all of them with the same period. Such a set is called
_isochronous_.

Let $x(t,z,\varepsilon)$ be the solution of system (2.2) such that
$x(0,z,\varepsilon)=z$. We write the linearization of the unperturbed system
(2.1) along the solution $x(t,z,0)$ as

$\dot{y}=D_{x}F_{0}(t,x(t,z,0))y,$ (2.3)

and we denote by $\mathcal{M}_{z}(t)$ the fundamental matrix of the linear
differential system (2.3) such that $\mathcal{M}_{z}(0)$ is the $n\times n$
identity matrix, and by $\xi:M=M_{k}\times M_{n-k}\to M_{k}$ the projection of
$M$ onto its first $k$ coordinates, that is,
$\xi(x_{1},\dots,x_{n})=(x_{1},\dots,x_{k})$.

The following results give sufficient conditions for the existence of limit
cycles for a system of the form (2.2) bifurcating from the periodic orbits of
system (2.1).

######  Theorem 4.

Let $V\subseteq M_{k}$ be an open and bounded set, and let
$\beta_{0}:\overline{V}\to M_{n-k}$ be a $C^{2}$ function. Assume

(i)
$\mathcal{Z}=\\{z_{\alpha}=(\alpha,\beta_{0}(\alpha))\,:\,\alpha\in\overline{V}\\}\subset
M$ and for each $z_{\alpha}\in\mathcal{Z}$ the solution $x(t,z_{\alpha},0)$ of
system (2.1) is T-periodic.

(ii) For each $z_{\alpha}\in\mathcal{Z}$, there is a fundamental matrix
$\mathcal{M}_{z_{\alpha}}(t)$ of system (2.3) such that the matrix
$\mathcal{M}_{z_{\alpha}}^{-1}(0)-\mathcal{M}_{z_{\alpha}}^{-1}(T)$ has the
$k\times(n-k)$ zero matrix in the upper right corner, and a $(n-k)\times(n-k)$
matrix $\Delta_{\alpha}$ in the lower right corner with
$\det(\Delta_{\alpha})\neq 0$.

Consider the function $\mathcal{F}:\overline{V}\to\mathbb{R}^{k}$ defined by

$\mathcal{F}(\alpha)=\xi\bigg{(}\int_{0}^{T}\mathcal{M}_{z_{\alpha}}^{-1}(t)F_{1}(t,x(t,z_{\alpha},0))\,dt\bigg{)}.$

If there exists $a\in V$ with $\mathcal{F}(a)=0$ and with
$\det(\mathcal{DF}(a))\neq 0$, then there is a limit cycle $x(t,\varepsilon)$
of period $T$ of system (2.2) such that $x(0,\varepsilon)\to z_{a}$ as
$\varepsilon\to 0$.

The result given by Theorem 4 can be found in the books of Malkin [15] and
Rosseau [16]. For a shorter proof, see [2]. There the result is proved in
$\mathbb{R}^{n}$, but it can be easily extended to a manifold $M$.

The next result allows to determine the existence of limit cycles in a system
of the form (2.2) in the case when there exists an open set, $V\subset M$,
such that for all $z\in\overline{V}$, the solution $x(t,z,0)$ is $T$-periodic.

######  Theorem 5.

Let $V\subseteq M$ be an open and bounded set with $\overline{V}\subseteq M$,
and assume that for all $z\in\overline{V}$ the solution $x(t,z,0)$ of system
(2.2) is $T$-periodic. Consider the function
$\mathcal{F}:\overline{V}\to\mathbb{R}^{n}$ defined by

$\mathcal{F}(z)=\int_{0}^{T}M_{z}^{-1}(t)F_{1}(t,x(t,z,0))\,dt.$

If there exists $a\in V$ with $\mathcal{F}(a)=0$ and with
$\det(\mathcal{DF}(a))\neq 0$, then there is a limit cycle $x(t,\varepsilon)$
of period $T$ of system (2.2) such that $x(0,\varepsilon)\to a$ as
$\varepsilon\to 0$.

For the proof of Theorem 5 see Corollary $1$ of [2].

##  3 Proof of Theorem 1

We use the result from averaging theory given in Theorem 4 to deduce the
existence of a limit cycle of system (LABEL:pertorb-lineal-2), for some
$\varepsilon>0$ small enough, bifurcating from a periodic orbit of the same
system with $\varepsilon=0$.

Since the general solution of the differential system (1.1), corresponding to
system (LABEL:pertorb-lineal-2) with $\varepsilon=0$, is given by

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad
r(t)=(r_{0}-1)e^{t}+1,$

it is clear that all the periodic solutions of that system are parameterized
by

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad r(t)=1,$

with
$(\theta_{0},\varphi_{0})\in\mathbb{S}^{2}\backslash\\{\varphi_{0}=\pm\pi/2\\}$.
Then, all the periodic solutions have period $2\pi$ and they fill the
invariant sphere $\\{r=1\\}$ except for the poles, which are equilibrium
points.

Therefore, for applying Theorem 4 we take $M=\mathbb{S}^{2}\times\mathbb{R}$
and

$\displaystyle\begin{split}k&=2,\,n=3,\\\
M_{k}&=M_{2}=\\{(\theta,\varphi,r)\in M\,:\,r=1\\}\cong\mathbb{S}^{2},\\\
x&=(\theta,\varphi,r),\\\ \alpha&=(\theta_{0},\varphi_{0}),\\\
\beta_{0}(\alpha)&=\beta_{0}(\theta_{0},\varphi_{0})=1,\\\
z_{\alpha}&=(\alpha,\beta_{0}(\alpha))=(\theta_{0},\varphi_{0},1),\\\
V&=\\{(\theta,\varphi,r)\in
M\,:\,r=1,\,\varphi\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0})\\}\\\
&\quad\text{ with $\delta_{0}>0$ small enough such that }\\\
&\quad\,\varphi^{*}:=\frac{b_{1}(a_{0}+a_{3}+a_{1}\pi)-a_{1}(b_{0}+b_{3}+b_{2}\pi)}{a_{1}b_{2}-a_{2}b_{1}}\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0}),\\\
\mathcal{Z}&=\overline{V}\times\\{r=1\\},\\\
x(t,z_{\alpha},0)&=(\theta_{0}+t,\varphi_{0},1),\\\ F_{0}(t,x)&=(1,0,r-1),\\\
F_{1}(t,x)&=(a_{0}+a_{1}\theta+a_{2}\varphi+a_{3}r,\,b_{0}+b_{1}\theta+b_{2}\varphi+b_{3}r,\,c_{0}+c_{1}\theta+c_{2}\varphi+c_{3}r),\\\
F_{2}(t,x,\varepsilon)&=0,\\\ T&=2\pi,\end{split}$ (3.1)

where we took $V\subset M_{2}$ as an open subset that contains the periodic
orbit for which it bifurcates a limit cycle, as we shall see next.

The fundamental matrix $\mathcal{M}_{z_{\alpha}}(t)$ with
$\mathcal{M}_{z_{\alpha}}(0)=Id$ of system (2.3) with $F_{0}$ and
$x(t,z_{\alpha},0)$ described above is the matrix
$\mathcal{M}_{z_{\alpha}}(t)=\exp(D_{x}F_{0}\,t)$, _i.e._

$\mathcal{M}_{z_{\alpha}}(t)=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&e^{t}\end{pmatrix}.$

Note that since $F_{0}$ defines a linear differential system, the fundamental
matrix $\mathcal{M}_{z_{\alpha}}(t)$ is independent of the initial conditions
$z_{\alpha}$.

We also have

$\mathcal{M}_{z_{\alpha}}^{-1}(0)-\mathcal{M}_{z_{\alpha}}^{-1}(2\pi)=\begin{pmatrix}0&0&0\\\
0&0&0\\\ 0&0&1-e^{-2\pi}\end{pmatrix},$

and therefore, all the assumptions in the in the statement of Theorem 4 are
satisfied.

With the described setting, the function
$\mathcal{F}(\alpha)=\mathcal{F}(\theta_{0},\varphi_{0})$ from the statement
of Theorem 4 associated with system (LABEL:pertorb-lineal-2) is

$\displaystyle\mathcal{F}(\theta_{0},\varphi_{0})$
$\displaystyle=\xi\bigg{(}\int_{0}^{2\pi}\mathcal{M}_{z_{\alpha}}^{-1}(t)F_{1}(\theta_{0}+t,\varphi_{0},1)\,dt\bigg{)}$
$\displaystyle=2\pi(a_{0}+a_{1}(\theta_{0}+\pi)+a_{2}\varphi_{0}+b_{3},\,b_{0}+b_{1}(\theta_{0}+\pi)+b_{2}\varphi_{0}+b_{3}).$

We have $\det(D\mathcal{F})=4\pi^{2}(a_{1}b_{2}-a_{2}b_{1})$, and therefore
$\det(D\mathcal{F})\neq 0$ for all $(\theta_{0},\varphi_{0})\in V$. Thus, the
only solution of $\mathcal{F}=0$ is given by

$\displaystyle\begin{split}\theta_{0}&=\frac{a_{2}(b_{0}+b_{3}+b_{1}\pi)-b_{2}(a_{0}+a_{3}+a_{1}\pi)}{a_{1}b_{2}-a_{2}b_{1}},\\\
\varphi_{0}&=\frac{b_{1}(a_{0}+a_{3}+a_{1}\pi)-a_{1}(b_{0}+b_{3}+b_{2}\pi)}{a_{1}b_{2}-a_{2}b_{1}}.\end{split}$
(3.2)

Note that such solution $(\theta_{0},\varphi_{0})$, where
$\varphi_{0}=\varphi^{*}$, is contained in the set $V$ described in (3.1).

Hence, by Theorem 4, if $\varepsilon>0$ is small enough, there is a periodic
solution, $(\theta(t,\varepsilon),\varphi(t,\varepsilon),r(t,\varepsilon))$,
of system (1.3), which is a limit cycle, and such that

$(\theta(0,\varepsilon),\varphi(0,\varepsilon),r(0,\varepsilon))\to(\theta_{0},\varphi_{0},1),$

when $\varepsilon\to 0$, and where $\theta_{0}$ and $\varphi_{0}$ are given in
(3.2).

##  4 Proof of Theorem 2

We use the result from averaging theory given in Theorem 4 to prove that, for
some $\varepsilon>0$ small enough, there exist a limit cycle of system
(LABEL:pertorb-lineal-high) bifurcating from a periodic orbit of the same
system with $\varepsilon=0$.

Since the general solution of system (LABEL:pertorb-lineal-high) with
$\varepsilon=0$ (that is, the one of system (1.4)), is

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad\nu(t)=\nu_{0}+t,\quad\phi(t)=\phi_{0},\quad
r(t)=(r_{0}-1)e^{t}+1,$

then all the periodic solutions of that system are

$\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},\quad\nu(t)=\nu_{0}+t,\quad\phi(t)=\phi_{0},\quad
r(t)=1,$

with
$(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})\in\mathbb{S}^{2}\backslash\\{\varphi_{0}=\pm\pi/2\\}\times\mathbb{S}^{2}\backslash\\{\varphi_{0}=\pm\pi/2\\}$.
That is, the periodic solutions fill the invariant manifold $\\{r=1\\}$ except
for the four equilibrium points $\\{\varphi=\pm\pi/2,\,\phi=\pm\pi/2\\}$, and
they have all period $2\pi$.

For applying Theorem 4 we take $M=(\mathbb{S}^{2})^{2}\times\mathbb{R}$ and

$\displaystyle\begin{split}k&=4,\,n=5,\\\
M_{k}&=M_{4}=\\{\theta,\varphi,\nu,\phi,r\in
M\,:\,r=1\\}\cong(\mathbb{S}^{2})^{2},\\\ x&=(\theta,\varphi,\nu,\phi,r),\\\
\alpha&=(\theta_{0},\varphi_{0},\nu_{0},\phi_{0}),\\\
\beta_{0}(\alpha)&=\beta_{0}(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})=1,\\\
z_{\alpha}&=(\alpha,\beta_{0}(\alpha))=(\theta_{0},\varphi_{0},\nu_{0},\phi_{0},1),\\\
V&=\\{(\theta,\varphi,\nu,\phi,r)\in
M\,:\,r=1,\,\varphi\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0})\\}\\\
&\quad\text{ with $\delta_{0}>0$ small enough such that
$\varphi_{0},\,\phi_{0}$ satisfying \eqref{sist-lin-cycle} satisfy}\\\
&\quad\,\varphi_{0},\,\phi_{0}\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0}),\\\
\mathcal{Z}&=\overline{V}\times\\{r=1\\},\\\
x(t,z_{\alpha},0)&=(\theta_{0}+t,\varphi_{0},\nu_{0}+t,\phi_{0},1),\\\
F_{0}(t,x)&=(1,0,1,0,r-1),\\\
F_{1}(t,x)&=\begin{pmatrix}a_{0}+a_{1}\theta+a_{2}\varphi+a_{3}\nu+a_{4}\phi+a_{5}r\\\
b_{0}+b_{1}\theta+b_{2}\varphi+b_{3}\nu+b_{4}\phi+b_{5}r\\\
c_{0}+c_{1}\theta+c_{2}\varphi+c_{3}\nu+c_{4}\phi+c_{5}r\\\
d_{0}+d_{1}\theta+d_{2}\varphi+d_{3}\nu+d_{4}\phi+d_{5}r\\\
e_{0}+e_{1}\theta+e_{2}\varphi+e_{3}\nu+e_{4}\phi+e_{5}r\end{pmatrix},\\\
F_{2}(t,x,\varepsilon)&=0,\\\ T&=2\pi,\end{split}$ (4.1)

where we chose $V\subset M_{4}$ as an open subset that contains the periodic
orbit for which it bifurcates a limit cycle, as we shall see next.

The fundamental matrix $\mathcal{M}_{z_{\alpha}}(t)$ with
$\mathcal{M}_{z_{\alpha}}(0)=Id$, of system (2.3) with $F_{0}$ and
$x(t,z_{\alpha},0)$ described above is the matrix
$\mathcal{M}_{z_{\alpha}}(t)=\exp(D_{x}F_{0}\,t)$, _i.e._ ,

$\mathcal{M}_{z_{\alpha}}(t)=\begin{pmatrix}1&0&0&0&0\\\ 0&1&0&0&0\\\
0&0&1&0&0\\\ 0&0&0&1&0\\\ 0&0&0&0&e^{t}\end{pmatrix}.$

We also have

$\mathcal{M}_{z_{\alpha}}^{-1}(0)-\mathcal{M}_{z_{\alpha}}^{-1}(2\pi)=\begin{pmatrix}0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&1-e^{-2\pi}\end{pmatrix},$

and therefore, all the assumptions in the statement of Theorem 4 are
satisfied.

With the described setting, the function
$\mathcal{F}(\alpha)=\mathcal{F}(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})$ in
the statement of Theorem 4 associated with system (LABEL:pertorb-lineal-high)
is

$\displaystyle\mathcal{F}(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})$
$\displaystyle=\xi\bigg{(}\int_{0}^{2\pi}\mathcal{M}_{z_{\alpha}}^{-1}(t)F_{1}(\theta_{0}+t,\varphi_{0},\nu_{0}+t,\phi_{0},1)\,dt\bigg{)}=(\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3},\mathcal{F}_{4}),$

with

$\displaystyle\mathcal{F}_{1}$
$\displaystyle=2\pi(a_{0}+a_{1}\theta_{0}+a_{1}\pi+a_{2}\varphi_{0}+a_{3}\nu_{0}+a_{3}\pi+a_{4}\phi_{0}+a_{5}),$
$\displaystyle\mathcal{F}_{2}$
$\displaystyle=2\pi(b_{0}+b_{1}\theta_{0}+b_{1}\pi+b_{2}\varphi_{0}+b_{3}\nu_{0}+b_{3}\pi+b_{4}\phi_{0}+b_{5}),$
$\displaystyle\mathcal{F}_{3}$
$\displaystyle=2\pi(c_{0}+c_{1}\theta_{0}+c_{1}\pi+c_{2}\varphi_{0}+c_{3}\nu_{0}+c_{3}\pi+c_{4}\phi_{0}+c_{5}),$
$\displaystyle\mathcal{F}_{4}$
$\displaystyle=2\pi(d_{0}+d_{1}\theta_{0}+d_{1}\pi+d_{2}\varphi_{0}+d_{3}\nu_{0}+d_{3}\pi+d_{4}\phi_{0}+d_{5}).$

Also, we have

$\det(D\mathcal{F})=16\pi^{4}\,\det\begin{pmatrix}a_{1}&a_{2}&a_{3}&a_{4}\\\
b_{1}&b_{2}&b_{3}&b_{4}\\\ c_{1}&c_{2}&c_{3}&c_{4}\\\
d_{1}&d_{2}&d_{3}&d_{4}\end{pmatrix}\neq 0,$

by assumption. The initial conditions
$(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})$ such that
$\mathcal{F}(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})=0$ are the solutions of
the linear system

$\displaystyle\begin{split}&a_{1}\theta_{0}+a_{2}\varphi_{0}+a_{3}\nu_{0}+a_{4}\phi_{0}=-a_{0}-a_{1}\pi-
a_{3}\pi-a_{5},\\\
&b_{1}\theta_{0}+b_{2}\varphi_{0}+b_{3}\nu_{0}+b_{4}\phi_{0}=-b_{0}-b_{1}\pi-
b_{3}\pi-b_{5},\\\
&c_{1}\theta_{0}+c_{2}\varphi_{0}+c_{3}\nu_{0}+c_{4}\phi_{0}=-c_{0}-c_{1}\pi-
c_{3}\pi-c_{5},\\\
&d_{1}\theta_{0}+d_{2}\varphi_{0}+d_{3}\nu_{0}+d_{4}\phi_{0}=-d_{0}-d_{1}\pi-
d_{3}\pi-d_{5}.\\\ \end{split}$ (4.2)

Since $\det(D\mathcal{F})\neq 0$, system (LABEL:sist-lin-cycle) has a unique
solution. Note that such solution $(\theta_{0},\varphi_{0},\nu_{0},\phi_{0})$
is contained in the set $V$ described in (4.1).

Hence, by Theorem 4, if $\varepsilon>0$ is small enough, there is a periodic
solution,

$(\theta(t,\varepsilon),\varphi(t,\varepsilon),\nu(t,\varepsilon),\phi(t,\varepsilon),r(t,\varepsilon)),$

of system (LABEL:pertorb-lineal-high), which is a limit cycle, and such that

$(\theta(0,\varepsilon),\varphi(0,\varepsilon),\nu(0,\varepsilon),\phi(0,\varepsilon),r(0,\varepsilon))\to(\theta_{0},\varphi_{0},\nu_{0},\phi_{0},1),$

when $\varepsilon\to 0$, and where $\theta_{0},\,\varphi_{0},\,\nu_{0},$ and
$\phi_{0}$ are given by the unique solution of system (LABEL:sist-lin-cycle).

##  5 Proof of Theorem 3

Since the general solution of system (LABEL:pertorb-lineal-r2s2) with
$\varepsilon=0$ is given by

$x(t)=x_{0}\cos(t)-y_{0}\sin(t),\quad
y(t)=x_{0}\sin(t)+y_{0}\cos(t),\quad\theta(t)=\theta_{0}+t,\quad\varphi(t)=\varphi_{0},$

the whole phase space is filled by periodic solutions, except form the
equilibrium points $(x,y,\theta,\varphi)=(0,0,\theta,-\pi/2)$ and
$(x,y,\theta,\varphi)=(0,0,\theta,\pi/2)$. Hence, the periodic solutions of
the differential system (1.6) fill an open set of the phase space
$M=\mathbb{R}^{2}\times\mathbb{S}^{2}$.

To prove Theorem 3 we use the result given in Theorem 5 to deduce that there
exist a limit cycle of system (LABEL:pertorb-lineal-r2s2), for some
$\varepsilon>0$ small enough, bifurcating from the periodic orbits of the same
system with $\varepsilon=0$.

To clarify the notation, here the solution $x(t,z,0)$ from the statement of
Theorem 5 will be denoted by $\mathbf{x}(t,z,0)$, and $x$ will denote the
first variable in the phase space.

For applying Theorem 5 we take $M=\mathbb{R}^{2}\times\mathbb{S}^{2}$ and

$\displaystyle\begin{split}\mathbf{x}&=(x,y,\theta,\varphi),\\\
z&=(x_{0},y_{0},\theta_{0},\varphi_{0}),\\\
\mathbf{x}(t,z,0)&=(x(t),y(t),\theta(t),\varphi(t))\text{ given by \eqref{sol-
th-r2s2}}\\\ F_{0}(t,x)&=(-y,x,1,0),\\\
F_{1}(t,x)&=\begin{pmatrix}a_{0}+a_{1}x+a_{2}y+a_{3}\theta+a_{4}\varphi\\\
b_{0}+b_{1}x+b_{2}y+b_{3}\theta+b_{4}\varphi\\\
c_{0}+c_{1}x+c_{2}y+c_{3}\theta+c_{4}\varphi\\\
d_{0}+d_{1}x+d_{2}y+d_{3}\theta+d_{4}\varphi\\\ \end{pmatrix},\\\
F_{2}(t,x,\varepsilon)&=0,\\\ T&=2\pi,\\\ V&=\\{(x,y,\theta,\varphi)\in
M\,:\,\|(x,y)\|<1+\kappa,\,\varphi\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0})\\},\\\
&\quad\ \text{with }\
\kappa=\frac{2\,\sqrt{a_{3}^{2}+b_{3}^{2}}}{\sqrt{a_{1}^{2}+a_{2}^{2}+b_{1}^{2}+b_{2}^{2}+2a_{1}b_{2}-2a_{2}b_{1}}},\\\
&\quad\text{ and with $\delta_{0}>0$ small enough such that }\\\
&\quad\,\varphi^{*}:=\frac{c_{0}d_{3}-c_{3}d_{0}}{c_{3}d_{4}-c_{4}d_{3}}\in(-\tfrac{\pi}{2}+\delta_{0},\tfrac{\pi}{2}-\delta_{0}).\end{split}$
(5.1)

where we chose $V\subset M$ as an open subset that contains the periodic orbit
for which it bifurcates a limit cycle, as we shall see next.

The fundamental matrix $\mathcal{M}_{z}(t)$ of system (2.3) with
$\mathcal{M}_{z}(0)=Id$ and with $F_{0}$ and $\mathbf{x}(t,z,0)$ described in
(5.1) is given by

$M_{z}(t)=\begin{pmatrix}\cos(t)&-\sin(t)&0&0\\\ \sin(t)&\cos(t)&0&0\\\
0&0&1&0\\\ 0&0&0&1\end{pmatrix}.$

Therefore all the assumptions in the statement of Theorem 5 are satisfied.

With the described setting the function
$\mathcal{F}(z)=\mathcal{F}(x_{0},y_{0},\theta_{0},\varphi_{0})$ in the
statement of Theorem 5 associated with system (LABEL:pertorb-lineal-r2s2),
namely,

$\displaystyle\mathcal{F}(x_{0},y_{0},\theta_{0},\varphi_{0})$
$\displaystyle=\int_{0}^{2\pi}M_{z}^{-1}(t)F_{1}(t,\mathbf{x}(t,x,0))\,dt,$

is given by
$\mathcal{F}=(\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3},\mathcal{F}_{4})$,
which after some straightforward computations can be written as

$\displaystyle\mathcal{F}_{1}$ $\displaystyle=(\pi a_{2}-\pi b_{1})y_{0}+(\pi
b_{2}+\pi a_{1})x_{0}-2\pi b_{3},$ $\displaystyle\mathcal{F}_{2}$
$\displaystyle=(\pi b_{2}+\pi a_{1})y_{0}+(\pi b_{1}-\pi a_{2})x_{0}+2\pi
a_{3},$ $\displaystyle\mathcal{F}_{3}$ $\displaystyle=2\pi
c_{3}\theta_{0}+2\pi c_{4}\varphi_{0}+2\pi^{2}c_{3}+2\pi c_{0},$
$\displaystyle\mathcal{F}_{4}$ $\displaystyle=2\pi d_{3}\theta_{0}+2\pi
d_{4}\varphi_{0}+2\pi^{2}d_{3}+2\pi d_{0}.$

Assuming that

$\det(\mathcal{DF})=\det\begin{pmatrix}\pi(b_{2}+a_{1})&\pi(a_{2}-b_{1})&0&0\\\
\pi(b_{1}-a_{2})&\pi(b_{2}+a_{1})&0&0\\\ 0&0&2\pi c_{3}&2\pi c_{4}\\\ 0&0&2\pi
d_{3}&2\pi d_{4}\\\ \end{pmatrix}\neq 0,$ (5.2)

the linear system
$(\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3},\mathcal{F}_{4})=(0,0,0,0)$
has a unique solution, given by

$\displaystyle\begin{split}x_{0}&=\frac{(2b_{2}+2a_{1})b_{3}-2a_{3}b_{1}+2a_{2}a_{3}}{b_{2}^{2}+b_{1}^{2}+a_{2}^{2}+a_{1}^{2}+2a_{1}b_{2}-2a_{2}b_{1}},\\\
y_{0}&=-\frac{(2b_{1}-2a_{2})b_{3}+2a_{3}b_{2}+2a_{1}a_{3}}{b_{2}^{2}+b_{1}^{2}+a_{2}^{2}+a_{1}^{2}+2a_{1}b_{2}-2a_{2}b_{1}},\\\
\theta_{0}&=-\frac{(\pi c_{3}+c_{0})d_{4}-\pi
c_{4}d_{3}-c_{4}d_{0}}{c_{3}d_{4}-c_{4}d_{3}},\\\
\varphi_{0}&=\frac{c_{0}d_{3}-c_{3}d_{0}}{c_{3}d_{4}-c_{4}d_{3}}.\end{split}$
(5.3)

Note that such solution $(x_{0},y_{0},\theta_{0},\varphi_{0})$, where
$\varphi_{0}=\varphi^{*}$, is contained in the set $V$ described in (5.1).

The condition (5.2) is clearly satisfied for all
$(x_{0},y_{0},\theta_{0},\varphi_{0})\in V$ taking into account the
assumptions in the statement of Theorem 3.

Hence, by Theorem 5, there is a periodic solution
$(x(t,\varepsilon),y(t,\varepsilon),\theta(t,\varepsilon),\varphi(t,\varepsilon))$
of system (LABEL:pertorb-lineal-r2s2), which is a limit cycle, and such that

$(x(0,\varepsilon),y(0,\varepsilon),\theta(0,\varepsilon),\varphi(0,\varepsilon),r(0,\varepsilon))\to(x_{0},y_{0},\theta_{0},\varphi_{0})$

when $\varepsilon\to 0$, and where $x_{0},\,y_{0},\,\theta_{0}$ and
$\varphi_{0}$ are given in (5.3).

Acknowledgements. This work is supported by the Ministerio de Ciencia,
Innovación y Universidades, Agencia Estatal de Investigación grants
PID2019-104658GB-I00 and BES-2017-081570, and the H2020 European Research
Council grant MSCA-RISE-2017-777911.

## References

  * [1] V. I. Arnold, Ordinary Differential Equations, _2nd printing of the 1992 edn_ (2006). Berlin, Springer, Universitext.

  * [2] A. J. Buica, J. P. Françoise and J. Llibre, Periodic solutions of nonlinear periodic differential systems with a small parameter, _Commun. Pure Appl. Anal_ 6 (2007), 103-111.
  * [3] C. Christopher and N. G. Lloyd, Small-amplitude limit cycles in polynomial Liénard systems, _Nonlinear Differ. Equ. Appl_ 3 (1996), 183-190.
  * [4] A. Ferragut, J. Llibre and M. A. Teixeira, Hyperbolic periodic orbits from the bifurcation of a four-dimensional nonlinear center, _Int. J. Bifurcation Chaos Appl. Sci. Eng._ 17 (2007), 2159-2167.
  * [5] H. Giacomini, J. Llibre and M. Viano, On the nonexistence, existence and uniqueness of limit cycles, _Nonlinearity_ 9 (1996), 501-516.
  * [6] M. Han and L. Li, Lower bounds for the Hilbert number of polynomial systems, _J. Differential Equations_ 252 (2012), 3278-3304.
  * [7] P. Hartman, A lemma in the theory of structural stability of differential equations, _Proc. Am. Math. Soc_ 11 (1960), 610-620.
  * [8] Yu. Ilyashenko, Centennial history of Hilbert’s 16th problem, _Bull. Am. Math. Soc_ 39 (2002), 301-354 (new series).
  * [9] G. Laferriere, G. J. Pappas and S. Yovine, Symbolic reachability computation for families of linear vector fields, _J. Symb. Comput._ 32 (2001), 231-253.
  * [10] G. Laferriere, G. J. Pappas and S. Yovine, A new class of decidable hybrid systems, _Hybrid Systems: Computation and Control (LNCS vol 1569)_ (1999) Berlin, Springer, 137-151.
  * [11] J. Llibre, R. M. Martins and M. A. Teixeira, Periodic orbits, invariant tori and cylinders of Hamiltonian systems near integrable ones having a return map equal to the identity. _J. Math. Phys_ 51 (2010), 082704.
  * [12] J. Llibre and M. A. Teixeira, Limit cycles bifurcating from a $2$-dimensional isochronous cylinder, _Appl. Math. Lett._ 22 (2009), 321-341.
  * [13] J. Llibre, M. A. Teixeira and J. Torregrosa, Limit cycles bifurcating form a $k$-dimensional isochronous center contained in $\mathbb{R}^{n}$ with $k\leq n$, _Phys. Anal. Geom._ 10 (2007), 237-249.
  * [14] J. Llibre and X. Zhang, Limit cycles of linear vector fields on manifolds, _Nonlinearity_ 29 (2016), 3120-3131.
  * [15] I. G. Malkin, Some problems of the theory of nonlinear oscillations (Russian), _Gosudarstv. Izdat. Tehn. - Theor. Lit._ (1956)
  * [16] M. Rosseau, Vibrations non linéaires et théorie de la stabilité (French), _Springer Tracts in Natural Philosophy_ , vol 8 (1996) Berlin, Springer.
  * [17] J. A. Sanders, F. Verhulst and J. Murdock, Averaging Methods in Nonlinear Dynamical Systems _Applied Mathematical Sciences)_ vol 59 2nd edn (2007) New York, Springer.
  * [18] J. Sotomayor, Liçoes de equaçoes diferenciais ordinárias,  _Euclid Project vol 11_ (1979) Rio de Janeiro, Instituto de Matemática Pura e Aplicada.
  * [19] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems (1996) Berlin, Springer, Universitext.