\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2018 (2018), No. 117, pp. 1--16.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2018 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2018/117\hfil Isochronous families of limit cycles]
{Isochronous families of limit cycles}

\author[R. Cobiaga, W. Reartes \hfil EJDE-2018/117\hfilneg]
{Romina Cobiaga, Walter Reartes}

\address{Romina Cobiaga \newline
Universidad Nacional del Sur,
 Av. Alem 1253, 8000 Bah\'ia Blanca,
Buenos Aires, Argentina}
\email{romina.cobiaga@uns.edu.ar}

\address{Walter Reartes \newline
Universidad Nacional del Sur,
 Av. Alem 1253, 8000 Bah\'ia Blanca,
Buenos Aires, Argentina}
\email{walter.reartes@gmail.com}

\thanks{Submitted March 30, 2017. Published May 15, 2018.}
\subjclass[2010]{34A34, 34C25}
\keywords{Limit cycle; isochronous orbit; delay differential equation}

\begin{abstract}
 In this article we present a method for determining if the frequency
 of a family of periodic orbits remains constant when a parameter changes.
 Two-dimensional systems of ordinary and delayed differential equations
 are considered. Several examples are given.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

In dynamic systems with oscillations, it is of great interest to study 
the frequency of such oscillations. Clearly this knowledge will be 
of importance in the development of applications in engineering and other
branches of technology.

For two-dimensional systems, the so-called isochronous centers have been 
studied in \cite{chavarriga99, cima99, gine07, mardesic95, mardesic06, sabatini05}. 
These are centers for which all orbits have the same frequency. 
Related to isochronous centers are isochronous foci 
\cite{gine03, gine07, sabatini03}. Here there are no periodic orbits and 
the time of return to a transversal cross-section to the orbit in the 
focus is studied.

The period of the emergent orbits of certain bifurcations has been studied 
in \cite{gasull05}. This is a particular case of a one-parameter dependent 
family of cycles. Such families have been studied for example in 
\cite{duff53, perko90}. Limit cycles of particular equations have also been 
studied, see  \cite{andersen82}.

In \cite{bel14} it has been proved that, for certain second order equations
 with delay, bifurcation branches appearing in a neighborhood of Hopf 
bifurcation due to the delay are isochronous.

A central problem is to know when the frequency of the oscillation does not 
depend on a parameter of the system. In this work, a methodology is 
developed to determine if a periodic orbit in two-dimensional systems of 
differential equations and differential equations with delay changes its 
frequency by varying a parameter. We call isochronous orbit with respect 
to the parameter an orbit that behaves this way. We also call isochronous 
family of limit cycles a family of this type.

In a similar spirit to that used in applying the theory of Floquet in 
determining the stability of a cycle, the problem of determining whether a
 periodic orbit is isochronous is reduced to studying a linear system with 
periodic coefficients. In this case, the existence of certain periodic solutions
of this new system, which we call auxiliary system, is equivalent to the isochronous 
character of the orbit.

\section{Auxiliary system}\label{cdb}

We begin by studying the case in which the periodic orbit is a limit cycle 
in a system of differential equations in the plane. Later, we will show the
necessary modifications to treat the case of differential equations with delay.

Consider the following system of ordinary differential equations
\begin{equation} \label{primera}
\begin{aligned}
\bar{x}_1' &= f_1(\bar{x}_1,\bar{x}_2,\mu),\\
\bar{x}_2' &= f_2(\bar{x}_1,\bar{x}_2,\mu).
\end{aligned}
\end{equation}
Here, $\mu$ is a real parameter, $\bar{x}_1$ and $\bar{x}_2$ are functions
of the variable $\bar{t}$, and the real functions $f_1$ and $f_2$ are
sufficiently differentiable. Suppose that this system has a limit
cycle with  frequency $\omega(\mu)>0$. Calling
$\mathbf{\bar{x}}_\mu=(\bar{x}_{1\mu},\bar{x}_{2\mu})$ the cycle and
 $\mathbf{f}=(f_1,f_2)$ the field, we have
\begin{equation} \label{solper}
\mathbf{\bar{x}}_\mu' = \mathbf{f}(\mathbf{\bar{x}}_\mu,\mu),\quad\text{where}\quad
\mathbf{\bar{x}}_\mu(\bar{t}) = \mathbf{\bar{x}}_\mu(\bar{t}+2\pi/\omega(\mu)).
\end{equation}

After a change of the independent variable to $t=\omega(\mu)\bar{t}$,
 \eqref{solper} becomes
\begin{equation}\label{frec1}
\omega(\mu)\mathbf{x}_\mu' = \mathbf{f}(\mathbf{x}_\mu,\mu),
\end{equation}
where the $2\pi$-periodic function ${\bf{x}}_\mu=(x_{1\mu},x_{2\mu})$ 
is continuous in the new variable $t$.

If we call $\Omega$ the space of continuous functions from the unit circle 
$S^1$ to $\mathbb{R}^2$, then the orbit that is solution of system 
\eqref{frec1} provides a continuous application, $\mu\mapsto\mathbf{x}_\mu$, 
of some interval $I\subset\mathbb{R}$ into $\Omega$. 
Without loss of generality, we can assume that the orbit surrounds the 
origin and suppose that it intersects the positive $x_1$-axis transversely
 at a point $\hat{x}_1(\mu)$. An additional condition must be set to define 
the $\mathbf{x}_\mu$ family due to the $S^1$-invariance of the periodic orbits.
 This can be done by imposing the condition that the orbits cross the $x_1$-axis 
at time $t=0$ (of course for this construction we can take any straight line 
that passes through the origin). We then have
\begin{equation}\label{condinic}
\mathbf{x}_\mu(0)=(\hat{x}_1(\mu),0).
\end{equation}
We assume that the function $\hat{x}_1(\mu)$ is differentiable with respect 
to $\mu$. Thus, $\mathbf{x}_\mu$ is continuous in both $\mu$ and $t$.

Differentiating \eqref{frec1} with respect to $\mu$ and assuming that the 
frequency is constant; i.e., $\omega(\mu)=\omega_0$, and therefore 
$\partial\omega/\partial\mu=0$, we obtain the following system of non-autonomous 
linear equations
\begin{equation}\label{xis}
\omega_0{\boldsymbol{\xi}_\mu'} 
= \frac{\partial\mathbf{f}}{\partial \mathbf{x}_\mu}{\boldsymbol{\xi}_\mu} 
+ \frac{\partial\mathbf{f}}{\partial\mu},
\end{equation}
where 
\begin{equation}\label{funxis}
\boldsymbol{\xi}_\mu(t)=(\xi_{1\mu}(t),\xi_{2\mu}(t))
=\Big(\frac{\partial x_{1\mu}}{\partial\mu}(t),
 \frac{\partial x_{2\mu}}{\partial\mu}(t)\Big),
\end{equation}
and $\partial\mathbf{f}/\partial\mathbf{x}_\mu$ the Jacobian matrix 
of $\mathbf{f}$ with respect to $x_{1\mu}$ and $x_{2\mu}$. 
The coefficients are periodic functions of period $2\pi$ because they are
 evaluated on the periodic orbit $\mathbf{x}_\mu$. We call this system the 
auxiliary system.

The variation of the limit cycle with $\mu$ is such that at $t=0$ 
it moves horizontally with velocity
\begin{equation}\label{derivmu}
{\boldsymbol{\xi}_\mu}(0) 
= \Big(\frac{\partial \hat{x}_1(\mu)}{\partial\mu},0\Big).
\end{equation}
Here, ${\boldsymbol{\xi}_\mu}(0)$ is on the horizontal axis by the 
particular choice \eqref{condinic} we have taken.

Next theorem ensures that the problem of studying the isochronous character 
of the family of limit cycles of equation \eqref{primera} is transformed 
into finding periodic solutions, of period $2\pi$, of system \eqref{xis}.

\begin{theorem} \label{thm2.1}
Suppose that system \eqref{primera} has a limit cycle for all points of some
 interval $I\subset\mathbb{R}$ of the parameter $\mu$. Also, suppose that 
it intersects the horizontal axis transversely at $t=0$. Then, the family 
is isochronous with respect to $\mu$; i.e., the frequency does not depend 
on the parameter in the interval, if and only if for some $\mu=\mu_0$ 
the limit cycle has frequency $\omega_0$, and for all $\mu$ in $I$, 
system \eqref{xis} has a periodic solution of period $2\pi$, with initial 
condition \eqref{derivmu}.
\end{theorem}

\begin{proof}
Suppose that system \eqref{primera} has a family of isochronous cycles 
in some interval of the parameter $\mu$; that is, it has a periodic solution 
$\mathbf{\bar{x}_\mu}=(\bar{x}_{1\mu},\bar{x}_{2\mu})$ of constant 
frequency $\omega(\mu)=\omega_0$. Then, system \eqref{frec1} becomes
\begin{equation}\label{frec2}
\omega_0 \mathbf{x}_\mu' = \mathbf{f}(\mathbf{x}_\mu,\mu),
\end{equation}
with ${\bf{x}}_\mu=(x_{1\mu},x_{2\mu})$ periodic of period $2\pi$. 
Differentiating with respect to $\mu$, we obtain that functions 
\eqref{funxis} verify equations \eqref{xis}. 
Therefore, $\boldsymbol{\xi}_\mu=\partial \mathbf{x}_\mu/\partial\mu$ 
is the periodic solution of period $2\pi$ sought.

Conversely, suppose that, for each $\mu$, system \eqref{xis} has a solution 
${\boldsymbol{\xi}_\mu}$ of period $2\pi$ with initial condition \eqref{derivmu}.
 Take $\mu=\mu_0$ in the interval $I$ and $\mathbf{x}_{\mu_0}$ the periodic 
solution corresponding to that parameter value. Define $\mathbf{x}_\mu$ as follows:
\begin{equation}\label{intedexi}
\mathbf{x}_\mu = \mathbf{x}_{\mu_0}+\int_{\mu_0}^\mu {\boldsymbol{\xi}_\nu}\,d\nu.
\end{equation}
As can be easily verified, this $\mathbf{x}_\mu$ is solution of \eqref{frec1} 
with the appropriate initial conditions and with $\omega(\mu)=\omega_0$ constant.
\end{proof}

Once system \eqref{xis} is constructed, the next step is to find the periodic 
solution of period $2\pi$ that verifies condition \eqref{derivmu}. 
To find this orbit we proceed to compute a line determined by the end points,
 at time $t=2\pi$, of the trajectories beginning on the axis $\xi_1$. 
This line can provide us with a method to ensure that the desired limit 
cycle does not exist. If the  cycle exists, this line helps us to find 
it for a fixed value of $\mu$.

Linear system \eqref{xis} can be written as
\begin{equation}\label{sisxis}
\boldsymbol{\xi}'(t)=A(t)\boldsymbol{\xi}(t)+\mathbf{b}(t).
\end{equation}
If we know a principal fundamental matrix solution, $\Phi(t)$, 
of the associated homogeneous system we have the following expression 
for the solution
\begin{equation}\label{solxis}
\boldsymbol{\xi}(t)=\Phi(t)\boldsymbol{c}+\Phi(t)
\int_{0}^{t}{\Phi^{-1}(s)\mathbf{b}(s)}ds,
\end{equation}
where $\boldsymbol{c}$ is a constant initial condition.

As a consequence of linearity, the orbits that begin at the axis $\xi_1$; that is,
\begin{equation}\label{cixis}
\boldsymbol{\xi}(0)=\Phi(0)\boldsymbol{c}=\boldsymbol{c}=(\xi_{01},0),
\end{equation}
are found in time $t=2\pi$ on a same line in the plane $\xi_1$-$\xi_2$. 
Calling $(\xi_1(2\pi),\xi_2(2\pi))=(\xi ^*_1,\xi^*_2)$ the end point of 
the path in the plane $\xi_1$-$\xi_2$, these points obey the parametric 
equation as a function of $\xi_{01}$,
\begin{equation}\label{recta2}
\begin{aligned}
\xi^*_1 &= \Phi_{11}\xi_{01}+\gamma_1,\\
\xi^*_2 &= \Phi_{21}\xi_{01}+\gamma_2,
\end{aligned}
\end{equation}
where
\begin{equation}
(\gamma_1,\gamma_2)^T=\Phi(2\pi)\int_{0}^{2\pi}{\Phi^{-1}(s)\mathbf{b}(s)}ds,
\end{equation}
and $\Phi_{ij}$ are the components of the matrix $\Phi(2\pi)$.

Because all the trajectories of system \eqref{xis} that start at the 
$\xi_1$-axis end at the same line at $t=2\pi$, then by plotting two of 
them we can determine this line. In cases where it does not cut the $\xi_1$-axis, 
we can ensure that the family of cycles is not isochronous. 
This situation can occur when the line is horizontal at a nonzero value
 of the variable $\xi_2$ (see Figure \ref{fig_armonico}). 
However, crossing the $\xi_1$-axis does not guarantee the existence of 
periodic orbits as shown in the following examples. In this case, 
it is enough to find the orbit that starts from that crossing point to 
check its periodicity. In this work, all the trajectories of the auxiliary 
system have been found by means of numerical simulations.

\begin{remark} \rm
We note that in finding the periodic orbit in the $\xi_1$-$\xi_2$ plane we 
not only prove that the family of orbits is isochronous, but we also find 
how the cycle changes in the direction of the $x_1$-axis as $\mu$ varies. 
As stated above, this axis could be changed by anyone through the origin, 
with a suitable synchronization of the family of cycles. Moreover, a 
direction that does not pass through the origin can be chosen by using 
as a line of synchronization a differentiable curve that passes through 
the origin and divides the plane into exactly two regions.
\end{remark}

\begin{remark} \rm
The concept of isochronous family of limit cycles could be generalized, 
in principle, for the case where the cycles are in a center. 
In this case we would call them isochronous family of periodic orbits in 
a center. We must explain that this concept has nothing to do with that 
of isochronous center, where the different orbits, for a same value of 
the parameter, have the same frequency. In addition, a family of the type 
studied in this work can be built in almost any center, with the only condition 
that the frequency of the orbits varies from orbit to orbit, as will be seen below. 
This construction could be of importance in applications in control theory.

Suppose that a system as \eqref{primera} has a center. We will also suppose 
that the center is structurally stable. Without loss of generality,
 we can assume that the orbit is around the origin. Then, we can suppose, 
as in the case of limit cycles, that at $t=0$ a particular periodic orbit 
crosses the $x_1$-axis at the point $\hat{x}_1$. The frequency of this orbit
 depends on the parameter $\mu$ and on the point $\hat{x}_1$.
 We call $\mathbf{\bar{x}}_{\mu,\hat{x}_1}$ that orbit. It satisfies the equation
\begin{equation}\label{solpercentro}
\mathbf{\bar{x}}_{\mu,\hat{x}_1}' 
= \mathbf{f}(\mathbf{\bar{x}}_{\mu,\hat{x}_1},\mu),\quad\text{where}\quad
\mathbf{\bar{x}}_{\mu,\hat{x}_1}(\bar{t}) 
= \mathbf{\bar{x}}_{\mu,\hat{x}_1}(\bar{t}+2\pi/\omega).
\end{equation}

Analogously to the case of limit cycles, if we call $\omega(\mu,\hat{x}_1)$ 
to the frequency, then we can write
\begin{equation}\label{frec1centro}
\omega(\mu,\hat{x}_1)\mathbf{x}_{\mu,\hat{x}_1}' 
= \mathbf{f}(\mathbf{x}_{\mu,\hat{x}_1},\mu).
\end{equation}

The frequency of this cycle depends on the value of $\mu=\mu_0$ and on 
the point $\hat{x}_1=x_{10}$ where it crosses the $x_1$-axis. 
If the frequency changes by varying $\hat{x}_1$; that is,
\begin{equation}
\frac{\partial\omega}{\partial \hat{x}_1}(\mu_0,x_{10})\ne 0,
\end{equation}
then, in a neighborhood of $\mu_0$, we can find a function $\hat{x}_1(\mu)$ 
satisfying
\begin{equation}
\omega(\mu,\hat{x}_1(\mu))=\omega_0.
\end{equation}
A family of isochronous periodic orbits can be defined by varying $\hat{x}_1$ 
with the function $\hat{x}_1(\mu)$. This family gives the following auxiliary
equation
\begin{equation}\label{xiscentro}
\omega_0{\boldsymbol{\xi}_\mu'} 
= \frac{\partial\mathbf{f}}{\partial \mathbf{x}_{\mu,\hat{x}_1}}
{\boldsymbol{\xi}_\mu} + \frac{\partial\mathbf{f}}{\partial\mu}.
\end{equation}
Where 
\begin{equation}\label{funxiscentro}
\boldsymbol{\xi}_\mu(t)=(\xi_{1\mu}(t),\xi_{2\mu}(t))
=\Big(\frac{\partial x_{1\mu,\hat{x}_1}}{\partial\mu}(t),
\frac{\partial x_{2\mu,\hat{x}_1}}{\partial\mu}(t)\Big).
\end{equation}
Two examples will suffice to clarify these points. First we consider the 
linear center in a harmonic oscillator of equation
\begin{equation}\label{ecuarmonico}
\begin{aligned}
x_1' &= x_2,\\
x_2' &= -\mu^2 x_1.
\end{aligned}
\end{equation}
The frequency depends on $\mu$ but not on the initial condition $x_0$, 
therefore there is no family of isochronous cycles. Equation \eqref{xiscentro} 
reduces to
\begin{equation}\label{ecuarmonicoxis}
\begin{aligned}
\xi_{1\mu}' &= \xi_{2\mu},\\
\xi_{2\mu}' &= -\mu^2 \xi_{1\mu}-2x_{1\mu,x_0}\mu=-\mu^2 \xi_{1\mu}-2x_0\cos(\mu t).
\end{aligned}
\end{equation}
In Figure \ref{fig_armonico} we see that system \eqref{ecuarmonicoxis} 
(with $\mu=2$)  has no periodic orbits, since the line is horizontal and 
does not cross the $\xi_1$-axis.
\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth]{fig1} % armonico.pdf
\caption{Orbits of system \eqref{ecuarmonicoxis} for the harmonic 
oscillator with horizontal line $\xi_2=-6.28$.} \label{fig_armonico}
\end{center}
\end{figure}

In contrast, if we consider the anharmonic oscillator of equation
\begin{equation}\label{ecuanarmonico}
\begin{aligned}
x_1' &= x_2,\\
x_2' &= -x_1-\mu x_1^3,
\end{aligned}
\end{equation}
equation \eqref{xiscentro} reduces to
\begin{equation}\label{ecuanarmonicoxis}
\begin{aligned}
\xi_{1\mu}'  &= \xi_{2\mu},\\
\xi_{2\mu}'  &= -(1+3\mu x_{1\mu,x_0}^2)\xi_{1\mu}-x_{1\mu,x_0}^3.
\end{aligned}
\end{equation}
In Figure \ref{fig_anarmonico} we can see the presence of a periodic orbit 
of period $2\pi$ starting from $\xi_1=-1$. This negative value shows 
that the radius of the cycle decreases when $\mu$ increases.
\end{remark}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=.5\textwidth]{fig2} % anarmonico.pdf
\end{center}
\caption{Orbits of system \eqref{ecuanarmonicoxis} for the 
anharmonic oscillator with line $\xi_2=-11.726-11.726\xi_1$ and 
$\mu=1$.} \label{fig_anarmonico}
\end{figure}


\subsection{Systems with angular velocity independent of the parameter}\label{scvac}

One type of systems having a family of trivially isochronous cycles is one 
in which the angular velocity is constant and independent of $\mu$. 
For example, system
\begin{equation}\label{sishopfdef}
\begin{aligned}
x_1' &= \mu x_1 - x_2 - x_1 (2 x_1^2 + x_2^2),\\
x_2' &= x_1 + \mu x_2 - x_2 (2 x_1^2 + x_2^2),
\end{aligned}
\end{equation}
transformed to polar coordinates gives
\begin{equation}\label{polares}
\begin{aligned}
r' &= \mu r-\Big(\frac{3}{2}+\cos(2\theta)\Big)r^3,\\
\theta' &= 1.
\end{aligned}
\end{equation}
In this case equation \eqref{xis} becomes
\begin{equation}\label{sishopfdefxis}
\begin{aligned}
\xi_{1\mu}' &=  (-6 x_{1\mu}^2 - x_{2\mu}^2 + \mu) \xi_{1\mu}
  - (1 + 2 x_{1\mu} x_{2\mu}) \xi_{2\mu} + x_{1\mu},\\
\xi_{2\mu}' &=  (1 - 4 x_{1\mu} x_{2\mu}) \xi_{1\mu} 
 - (2 x_{1\mu}^2 + 3 x_{2\mu}^2 + \mu) \xi_{2\mu} + x_{2\mu}.
\end{aligned}
\end{equation}
In Figure \ref{figura1a} we see the cycle to which the trajectories that 
start from the $\xi_1$-axis converge ($\mu=1$ was taken).
 As it can be seen, line \eqref{recta2} coincides with the $\xi_1$-axis. 
This occurs in all cases of systems with constant angular velocity.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig3} % hopfdef.pdf
\end{center}
\caption{Orbits of system \eqref{xis} for equation \eqref{sishopfdef} with 
$\omega_0=1$ and $\mu=1$.} \label{figura1a}
\end{figure}

\subsection{van der Pol and Rayleigh equations}\label{rayyvan}

In this section we will show two examples of well-known equations whose 
families of cycles are not isochronous. However, the frequency of these 
cycles is known as a function of the parameter with arbitrary precision. 
Therefore, a change of variables can be made to obtain systems where 
the limit cycles are isochronous.

As it is known, the limit cycle of van der Pol equation
\begin{equation}\label{sisvdp}
\begin{aligned}
x_1' = x_2,\\
x_2' = -\mu(x_1^2-1)x_2-x_1
\end{aligned}
\end{equation}
is not isochronous with respect to $\mu$. Here, equation \eqref{xis} 
takes the form
\begin{equation}\label{sisvdpxis}
\begin{aligned}
\omega_0\xi_{1\mu}' &= \xi_{2\mu},\\
\omega_0\xi_{2\mu}' &=  -(1 + 2\mu x_{1\mu} x_{2\mu}) \xi_{1\mu} + \mu (1 - x_{1\mu}^2) \xi_{2\mu} + (1 - x_{1\mu}^2) x_{2\mu}(t).
\end{aligned}
\end{equation}
Several orbits of system \eqref{sisvdpxis} are shown in Figure \ref{figura1b}. 
Observing its behavior and the line given by equation \eqref{recta2} 
we see that the system does not have periodic orbits of period $2\pi$. 
The green curve starts at the point where the line intersects the $\xi_1$-axis, 
at the value $0.0134$, and clearly it is not a closed orbit of period $2\pi$.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig4} % vdpcompleto.pdf
\end{center}
\caption{Orbits of system \eqref{sisvdpxis} for van der Pol equation with $\mu=1$.
 Straight line $\xi_2 = 3174.35 \xi_1-42.58$ is shown.} \label{figura1b}
\end{figure}

As stated before, the frequency $\omega(\mu)$ of the limit cycle is known 
with sufficient approximation in power series of the parameter 
$\mu$ \cite{andersen82}. Making a change on the variable $t$, we obtain
\begin{equation}\label{sisvdpomega}
\begin{aligned}
 x_1' &= \omega(\mu)^{-1}x_2,\\
 x_2' &= -\omega(\mu)^{-1}(\mu(x_1^2-1)x_2-x_1),
\end{aligned}
\end{equation}
where we take the approximation
\begin{equation}\label{frecvdp}
\omega(\mu)=1-\frac{1}{16}\mu^2 + \frac{17}{3072}\mu^4 
+ \frac{35}{884736}\mu^6 - \frac{678899}{5096079360}\mu^8 
+ \frac{28160413}{2293235712000}\mu^{10},
\end{equation}
and the names of the variables are preserved for simplicity. 
In this system the family of cycles is isochronous 
(within the approximation given by \eqref{frecvdp}). 
In Figure \ref{figura2b} we observe the limit cycle in the auxiliary 
system for $\mu=1$, it is a periodic orbit of period $2\pi$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=.5\textwidth]{fig5} % vdp_isocronico_detalle.pdf
\end{center}
\caption{Cycles of the auxiliary system for van der Pol equation
 \eqref{sisvdpomega}, with $\mu=1$.} \label{figura2b}
\end{figure}


Similarly, the limit cycles of the Rayleigh oscillator,
\begin{equation}\label{ecray}
\begin{aligned}
x_1' &= x_2,\\
x_2' &= -\mu (x_2^2-1) x_2 - x_1,
\end{aligned}
\end{equation}
are not isochronous either. In Figure \ref{figura1c} we can see the orbits
 of the auxiliary system for this case and the corresponding line given 
by equation \eqref{recta2}. This line intersects the $\xi_1$-axis in the 
value $0.1915$, but in the graph it is observed that the orbit starting 
from that point, the green curve, is not periodic with period $2\pi$.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig6} % rayleigh.pdf
\end{center}
\caption{Orbits of the auxiliary system for Rayleigh equation \eqref{ecray}
 with $\mu=1$. Straight line $\xi_2 = -112.5 + 587.59\xi_1$ is shown.} 
\label{figura1c}
\end{figure}


In the same way as with van der Pol equation, system \eqref{ecray} becomes
\begin{equation}\label{ecrayomega}
\begin{aligned}
x_1' &=\omega(\mu)^{-1} x_2,\\
x_2' &= -\omega(\mu)^{-1}(\mu (x_2^2-1) x_2 - x_1).
\end{aligned}
\end{equation}
The frequency of these limit cycles can be calculated by the method 
of Poincar\'e-Lindstedt \cite{verhulst90}. We have the fact that the 
formula for $\omega(\mu)$ is the same as \eqref{frecvdp}. 
Again the obtained system is isochronous. We see in Figure \ref{figura2a} 
the cycle of period $2\pi$ of the auxiliary system for $\mu=1$.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig7} % rayleigh_isocronico.pdf
\end{center}
\caption{Cycles of the auxiliary system for Rayleigh equation \eqref{ecrayomega}, 
with $\mu=1 $} \label{figura2a}
\end{figure}


\section{Delay differential equations}\label{retardo}

In \cite{bel14} the problem of isochronous cycles in delayed systems of 
a particular type was considered. The differential equations with delay 
are a special case of functional differential equations \cite{hale93}. 
The latter are of the form
\begin{equation}\label{hale}
\mathbf{x}'(t)=F(t,\mathbf{x}_t),
\end{equation}
where $\mathbf{x}(t)\in\mathbb{R}^n,$ $\mathbf{x}_t\in C([-\tau,0],\mathbb{R}^n)$ 
is the function $\mathbf{x}_t(\theta)=\mathbf{x}(t+\theta)$ and 
$F:D\subset\mathbb{R} \times C([-\tau,0],\mathbb{R}^n)\to \mathbb{R}^n$ 
is continuous.

In the case of differential equations with delay, the functional $F$ 
is of the form $F(t,\mathbf{x}_t)=h(t, \mathbf{x}(t), \mathbf{x}(t - \tau_1), 
\ldots, \mathbf{x}(t - \tau_k))$, where 
$h\colon U\subset\mathbb{R}\times\mathbb{R}^{nk}\to\mathbb{R}^n$ 
is continuous and there is a constant $\tau$ such that $0<\tau_i\le\tau$, 
$i=1,\cdots,k$.

In \cite{bel14} equations of form
\begin{equation}\label{sgret}
x'' + g(x,\beta) = \gamma f(x-x^\tau)
\end{equation}
are considered, where $x(t),\beta$ and $\gamma\in\mathbb{R}$, and
 $x^\tau(t)=x(t-\tau)$ (we use the superscript to avoid confusion with 
the function $\mathbf{x}_t$). Function $g$ is $C^1$ and $f$ is real analytic. 
In some cases a Hopf bifurcation occurs when any of the parameters $\gamma$ 
 or  $\beta$ changes. In the above work it was proved that the generic 
situation is that the frequency of the emergent limit cycles is constant 
in a neighborhood of the bifurcation.

In this work we consider differential equations with delay in the plane, 
of the type
\begin{equation}\label{scr}
\mathbf{{\bar{x}}}'=\mathbf{f}({\bf{\bar x}},\mathbf{\bar x}^\tau,\mu).
\end{equation}


As in Section \ref{cdb}, we assume that system has a limit cycle,
 $\mathbf{\bar x}_\mu$, with frequency $\omega(\mu)$. 
Then by suitable change of variables it can be written as
\begin{equation}\label{scrc2pciclo}
\omega(\mu)\mathbf{x}'_\mu=\mathbf{f}
(\mathbf{x}_\mu,\mathbf{x}_\mu^{\omega\tau},\mu).
\end{equation}

Differentiating with respect to $\mu$ and assuming that the frequency 
does not depend on $\mu$, we get the following delayed auxiliary system
\begin{equation}\label{xisdel}
\omega_0{\boldsymbol{\xi}_\mu'} = \frac{\partial\mathbf{f}}{\partial 
\mathbf{x}_\mu}{\boldsymbol{\xi}_\mu} 
+  \frac{\partial\mathbf{f}}{\partial \mathbf{x}^{\omega_0\tau}_\mu}
{\boldsymbol{\xi}^{\omega_0\tau}_\mu} +\frac{\partial\mathbf{f}}{\partial\mu},
\end{equation}
where the variations $\boldsymbol{\xi}_\mu$ are given again by 
equations \eqref{funxis} and $\partial\mathbf{f}/\partial\mathbf{x}_\mu$, 
$\partial\mathbf{f}/\partial\mathbf{x}^{\omega_0\tau}_\mu$ are the 
Jacobians with respect to the variables and retarded variables, respectively.

Of course the interpretation of solutions of system \eqref{xisdel} 
is much more involved than the case of an ordinary differential equation. 
The difference essentially lies in that the initial condition is in an 
infinite dimensional space and it is not a point in $\mathbb{R}^2$.
 We use a heuristic approach to study this problem.

Suppose that system \eqref{scr} has a unique periodic solution through 
the point $(\hat{x}_1,0)$. For example, from a Hopf bifurcation. 
Let us further suppose that with a small variation of $\mu$ another cycle 
is obtained. Suppose that these cycles are stable, then if we integrate 
equation \eqref{scrc2pciclo} with a constant function in the interval 
$[-\tau,0]$, at the point $(\hat{x}_1,0)$, then the orbit will probably 
converge to the new cycle. It is interesting to know whether, with 
this variation of $\mu$, the frequency of the new cycle remains constant.

System \eqref{xisdel} helps us to answer this question. If the points on 
the $\xi_1$-axis are used as the origin of orbits with a constant initial 
condition during the delay time, then these orbits can converge to periodic 
orbits of period $2\pi$. If this happens we may consider it as an indication 
that the system is isochronous. On the other hand, if it fails to find cycles 
of period $2\pi$ with constant initial conditions of the type here considered, 
it is a strong evidence that the family of limit cycles is not isochronous. 
The difficulty in interpreting the solutions of the auxiliary system in delayed 
equations is that it is impossible to cover all the initial conditions to test 
the existence of periodic orbits.

An interesting situation we have found is that there may be infinite orbits 
of period $2\pi$, as will be seen in the examples below. This happens as follows.
 Suppose that $\boldsymbol{\xi}_\mu$ is a solution of period $2\pi$ of 
the auxiliary system. Then the family $\boldsymbol{\xi}_\mu+\boldsymbol{\eta}_\mu$ 
will also be a solution of period $2\pi$ if $\boldsymbol{\eta}_\mu$ satisfies 
the homogeneous equation
\begin{equation}\label{etadel}
\omega_0{\boldsymbol{\eta}_\mu'} 
= \frac{\partial\mathbf{f}}{\partial \mathbf{x}_\mu}{\boldsymbol{\eta}_\mu} 
+ \frac{\partial\mathbf{f}}{\partial \mathbf{x}^{\omega_0\tau}_\mu}
{\boldsymbol{\eta}^{\omega_0\tau}_\mu}.
\end{equation}
If $\omega_0\tau$ is a multiple of $\pi$, then this equation has infinite 
proportional solutions of period $2\pi$. One of them is the one that 
allows us to write $\mathbf{x}_\mu$ as a function of $\mathbf{x}_{\mu_0}$ 
using equation \eqref{intedexi}; that is,
\[ % \label{intedexiret}
\mathbf{x}_\mu = \mathbf{x}_{\mu_0}+\int_{\mu_0}^\mu {\boldsymbol{\xi}_\nu}\,d\nu.
\]
We state the following theorem.

\begin{theorem} \label{thm3.1}
Suppose that system \eqref{scr} has a limit cycle for all points of some 
interval $I\subset\mathbb{R}$ of the parameter $\mu$. Also, suppose that 
it crosses the horizontal axis transversely at $t=0$. Then the family 
is isochronous with respect to $\mu$; i.e., the frequency does not 
depend on the parameter in the interval, if and only if for some
 $\mu=\mu_0$ the limit cycle has frequency $\omega_0$, and for all $\mu$ 
in $I$ system \eqref{xisdel} has a periodic solution of period $2\pi$
 which is equal to the derivative of $\mathbf{x}_\mu$ with respect to 
$\mu$ for $t\in[-\omega_0\tau,0]$.
\end{theorem}

\begin{remark} \rm
The presence of an infinite number of cycles of period $2\pi$ in the 
auxiliary system can also occur in systems of ordinary differential equations. 
Here $\boldsymbol{\eta}_\mu$ must verify equation
\begin{equation}\label{etasindel}
\omega_0{\boldsymbol{\eta}_\mu'}
 = \frac{\partial\mathbf{f}}{\partial \mathbf{x}_\mu}{\boldsymbol{\eta}_\mu}.
\end{equation}
According to Floquet's theory \cite{floquet1883} the fundamental matrix of 
solutions of this system is of the form $\boldsymbol{P}(t)e^{t\boldsymbol{R}}$,
 where $\boldsymbol{P}$ is a $2\pi$-periodic matrix and $\boldsymbol{R}$ 
is a constant matrix. That is, there may be periodic solutions of period $2\pi$ 
only if any of the eigenvalues of  $\boldsymbol{R}$ (characteristic values) 
is $2\pi i$.
\end{remark}

In what follows we will show three examples. In the first one, it is known 
that the limit cycles are not isochronous \cite{bel12a}. For the other two, 
it was proved in \cite{bel14} that the branches of Hopf bifurcation are isochronous.

\subsection{van der Pol with delayed feedback}\label{vdpcr}

We consider the equation studied in \cite{bel12a},
\begin{equation}\label{vdp_retardo}
x'' + (x^2-1)x' + x = \mu x^\tau.
\end{equation}
When the parameter $\mu$ is zero, we have a van der Pol equation.
 For $\mu$ sufficiently small, the limit cycle is maintained but it 
is deformed, generating a family which is not isochronous. 
These cycles of frequency $\omega(\mu)$ give rise to cycles of frequency $1$ 
of system
\begin{equation}\label{vdpretciclo}
\begin{aligned}
\omega(\mu)x_{1\mu}' &= x_{2\mu},\\
\omega(\mu)x_{2\mu}' &= -(x_{1\mu}^2-1)x_{2\mu} - x_{1\mu} 
+ \mu x_{1\mu}^{\omega(\mu)\tau}.
\end{aligned}
\end{equation}
The auxiliary system takes the form
\begin{equation}\label{vdpretxis}
\begin{aligned}
\omega_0\xi_{1\mu}' &= \xi_{2\mu},\\
\omega_0\xi_{2\mu}' &= -(1+2x_{1\mu}x_{2\mu})\xi_{1\mu} + (1-x_{1\mu}^2)\xi_{2\mu} 
+ \mu\xi_{1\mu}^{\omega_0\tau} + x_{1\mu}^{\omega_0\tau}.
\end{aligned}
\end{equation}
Hopf bifurcations also appear in this system by varying the parameter $\mu$. 
These branches are not isochronous either.

As can be seen in Figure \ref{fig_vdpret_orbita} the orbits of the auxiliary 
system \eqref{vdpretxis} are not bounded, because they move away from the origin.
 This is an indication of the non-isochronous character of the family. 
The orbit of the figure was calculated with an initial condition constant 
during the delay time at a point in the $\xi_1$-axis.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig8} % vdpretorbita.pdf
\end{center}
\caption{Orbit of the auxiliary system \eqref{vdpretxis} for $\mu=0.1$, 
$\tau=9$ and $\omega_0=0.99315$.} \label{fig_vdpret_orbita}
\end{figure}

\subsection{Rotating pendulum with delay}\label{retejem}

In \cite {bel14} the following differential equation with delay was studied:
\begin{equation}\label{pendulo_orig}
x'' + (\beta-\cos x)\sin x = \mu \sin(x-x^\tau).
\end{equation}
It represents the movement of a pendulum restricted to move on the surface 
of a vertical plane that rotates with constant angular velocity. 
A delayed feedback, which may be interpreted as a torque acting on the pivot, 
was added. This equation has two equilibria, one around $x=0$ and another 
at $x^*=\arccos\beta$. It was proved that the branches of Hopf bifurcations 
that appear due to the delay are isochronous, i.e., the emergent cycles of 
this bifurcation around the nonzero equilibrium have constant frequency. 
Hopf bifurcations occur at frequency $\omega=(2n+1)\pi/\tau$ for $n=0,1,\ldots$.

To apply the method developed in the present work we rewrite equation 
\eqref{pendulo_orig} around the non-zero equilibrium. 
After the change of independent variable, the $2\pi$-periodic orbits verify 
the equation
\begin{equation}\label{pendulo}
\begin{aligned}
\omega(\mu)x_{1\mu}' &= x_{2\mu},\\
\omega(\mu)x_{2\mu}' &= -(\beta - \cos(x_{1\mu}-x^*))\sin(x_{1\mu}-x^*) + \mu \sin(x_{1\mu}-x_{1\mu}^{\omega(\mu)\tau}).
\end{aligned}
\end{equation}

The auxiliary system for this problem is
\begin{equation}\label{penduloretxis}
\begin{aligned}
\omega_0\xi_{1\mu}' &= \xi_{2\mu},\\
\omega_0\xi_{2\mu}' 
&= (\cos2(x_{1\mu}-x^*) -\beta\cos(x_{1\mu}-x^*) 
 + \mu\cos(x_{1\mu}-x^*))\xi_{1\mu}\\
&\quad+ \mu\cos(x_{1\mu}-x^*)\xi_{1\mu}^{\omega_0\tau}.
\end{aligned}
\end{equation}

The orbits of system \eqref{penduloretxis} converge to cycles of 
period $2\pi$ for initial conditions constant during the delay time 
at a point in the $\xi_1$-axis. This provides a verification of the 
isochronous character of the Hopf cycles in this system. In this case 
the cycles are not unique, as explained in Section \ref{retardo}. 
As the delay is a multiple of $\pi$, infinite solutions of period $2\pi$ 
of equation \eqref{etadel} are obtained.

If we plot the cycles to which the solutions of \eqref{penduloretxis} 
converge and in a third axis the initial value of the coordinate $\xi_1$, 
where the integration starts, we obtain Figure \ref{figpenduloret}. 
The surface determined by all these cycles is a ruled surface as also 
shown in Figure \ref{figpenduloret}, where the lines are represented 
by points on the lines joining the corresponding cycles. 
It is a consequence of the linearity of the equations. 
If, instead of taking initial conditions during the delay time as mentioned above, 
other conditions are taken (for example a non-constant function of time 
or starting from another line), the same cycles are obtained.


\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.48\textwidth]{fig9a} % ciclos_pendulo_3d.pdf
 \includegraphics[width=0.48\textwidth]{fig9b} % pendulo_sup_reglada.pdf
\end{center}
\caption{Left: cycles of the auxiliary system \eqref{penduloretxis}
 for delayed rotatory pendulum, with $\mu=-1$, $\beta=0.5$, $\tau=2$ and 
$\omega_0=\pi/\tau$. Right: dotted right lines joining different cycles}
 \label{figpenduloret}
\end{figure}

\subsection{Anharmonic oscillator with delay}\label{aacr}

Another similar example is shown in \cite{bel14}. It is given by the 
following differential equation with delay
\begin{equation}\label{anar_orig}
x'' + x + \beta x^3 = \mu (x-x^\tau).
\end{equation}
As in the previous example this equation undergoes Hopf bifurcations 
when the parameter $\mu$ varies. The emergent cycles have frequency 
$\omega=(2n+1)\pi/\tau$ for $n = 0,1,\ldots$. 
After reparametrization, equation \eqref{anar_orig} gives
\begin{equation}\label{osciladorretardo}
\begin{aligned}
\omega(\mu)x_1' &= x_2,\quad\\
\omega(\mu)x_2' &= -\beta x_1^3 - x_1 + \mu (x_1-x_1^{\omega(\mu)\tau}),
\end{aligned}
\end{equation}
and the auxiliary equation is
\begin{equation}\label{anarmonretxis}
\begin{aligned}
\omega_0\xi_{1\mu}' &= \xi_{2\mu},\\
\omega_0\xi_{2\mu}' &= (\mu-3\beta x_{1\mu}^2-1)\xi_{1\mu} 
- \mu\xi_{1\mu}^{\omega_0\tau} + x_1-x_1^{\omega_0\tau}.
\end{aligned}
\end{equation}

The behavior of this system is very similar to the previous one. 
In Figure \ref{figanaret} we observe the cycles obtained by integrating 
equation \eqref{anarmonretxis} from different points on the $\xi_1$-axis 
(which is represented on the vertical axis). The lines connecting the 
cycles are also observed as dotted lines. Again this shows that the cycles 
in this system are isochronous.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.48\textwidth]{fig10a} % ciclos_anarmonico_3d.pdf
 \includegraphics[width=0.48\textwidth]{fig10b} % anarmonico_sup_reglada.pdf
\end{center}
\caption{Left: cycles of the auxiliary system \eqref{anarmonretxis} 
for delayed anharmonic oscillator, with $\mu=-0.4$, $\beta=-1$, 
$\tau=2.5$ and $\omega_0=\pi/\tau$. Right: dotted right lines joining 
different cycles.} \label{figanaret}
\end{figure}

\section{Conclusions}\label{conclu}

We have developed a methodology to study the frequency behavior of a family 
of limit cycles in differential equations and differential equations with 
delay in the plane. In particular, this methodology allows us to know if 
the family has constant frequency when the parameter changes. 
The problem is reduced to study the periodic solutions of non-homogeneous 
and non-autonomous linear differential equations.

While the method presented in this article works for systems in the plane, 
it could be extended to higher dimensions. However, its utility would be 
compromised by the fact that the search for periodic solutions of the 
auxiliary system should be done in one plane.

For the case of differential equations with delay the search for periodic 
cycles of the auxiliary equation is much more involved than in ordinary equations. 
In particular, the straight line through which the solutions pass after 
integration for time $2\pi$ corresponds to a linear subvariety in the Banach 
space of continuous functions and therefore it is not accessible with the 
graphic methods developed here. The problem would deserve further study, 
perhaps with the application of topological techniques such as degree 
theory \cite{lloyd78}.

\subsection*{Acknowledgments}
This work is supported by the Universidad Nacional del Sur (Grant no. PGI 24/L096).

\begin{thebibliography}{10}

\bibitem{andersen82} C.~M. Andersen, J.~F. Geer;
\newblock Power expansions for the frequency and period of the limit cycle of
  the van der pol equation.
\newblock {\em Siam J. Appl. Math.}, 42(3):678--693, 1982.

\bibitem{bel12a} A.~Bel, W.~Reartes;
\newblock The homotopy analysis method in bifurcation analysis of delay
  differential equations.
\newblock {\em International Journal of Bifurcation and Chaos}, 22(8), 2012.

\bibitem{bel14} A.~Bel, W.~Reartes;
\newblock Isochronous bifurcations in second-order delay differential
  equations.
\newblock {\em Electronic Journal of Differential Equations}, 2014(162):1--12,
  2014.

\bibitem{chavarriga99} J.~Chavarriga, M.~Sabatini;
\newblock A survey of isochronous centers.
\newblock {\em Qualitative Theory of Dynamical Systems}, 1:1--70, 1999.

\bibitem{cima99} A.~Cima, F.~Ma\~nosas, J.~Villadelprat;
\newblock Isochronicity for several classes of hamiltonian systems.
\newblock {\em Journal of Differential Equations}, 157:373--413, 1999.

\bibitem{duff53} G.~F.~D. Duff;
\newblock Limit-cycles and rotated vector fields.
\newblock {\em Annals of Mathematics}, 57(1):15--31, 1953.

\bibitem{floquet1883} G.~Floquet;
\newblock Sur les \'equations diff\'erentielles lin\'eaires \`a coefficients
  p\'eriodiques.
\newblock {\em Annales scientifiques de l'\'Ecole Normale Sup\'erieure},
  12:47--88, 1883.

\bibitem{gasull05} A.~Gasull, V.~Ma\~nosa,  J.~Villadelprat;
\newblock On the period of the limit cycles appearing in one-parameter
  bifurcations.
\newblock {\em J. Differential Equations}, 213:255--288, 2005.

\bibitem{gine03} J.~Gin\'e;
\newblock Isochronous foci for analytic differential systems.
\newblock {\em Internat. J. Bifur. Chaos}, 13(6):1617--1623, 2003.

\bibitem{gine07} J.~Gin\'e, M.~Grau;
\newblock Characterization of isochronous foci for planar analytic differential
  systems.
\newblock {\em Proceedings of the Royal Society of Edinburgh: Section A
  Mathematics}, 135(05):985--998, 2007.

\bibitem{hale93} J.~K. Hale, S.~M.~Verduyn Lunel;
\newblock {\em Introduction to Functional Differential Equations}, volume~99 of
  {\em Applied Mathematical Sciences}.
\newblock Springer--Verlag, 1993.

\bibitem{lloyd78} N.~G. Lloyd;
\newblock {\em Degree Theory}.
\newblock Cambridge University Press, 1978.

\bibitem{mardesic06} P.~Marde\v{s}i\'c, D.~Mar\'in, J.~Villadelprat;
\newblock The period function of reversible quadratic centers.
\newblock {\em J. Differential Equations}, 224:120--171, 2006.

\bibitem{mardesic95} P.~Marde\v{s}i\'c, P.~Rousseau, B.~Toni;
\newblock Linearization of isochronous centers.
\newblock {\em J. Differential Equations}, 121:67--108, 1995.

\bibitem{perko90} L.~M. Perko;
\newblock Global families of limit cycles of planar analytic systems.
\newblock {\em Transactions of the American Mathematical Society},
  322(2):627--656, 1990.

\bibitem{sabatini03} M.~Sabatini;
\newblock Non-periodic isochronous oscillations in plane differential systems.
\newblock {\em Ann. Mat. Pura Appl.}, 182(4):487--501, 2003.

\bibitem{sabatini05} M.~Sabatini;
\newblock On the period function of planar systems with unknown normalizers.
\newblock {\em Proceedings of the American Mathematical Society},
  134(2):531--539, 2005.

\bibitem{verhulst90} F.~Verhulst;
\newblock {\em Nonlinear Differential Equations and Dynamical Systems}.
\newblock Springer-Verlag, 1992.

\end{thebibliography}

\end{document}