text:
\begin{document}

\title{Distinguished trajectories in time dependent vector fields}

\author{J. A. Jim\'enez Madrid, Ana M. Mancho\\
Instituto de Ciencias Matem\'aticas, CSIC-UAM-UC3M-UCM, \\
Serrano 121, 28006 Madrid (Spain)}

\maketitle

\begin{abstract}

We introduce a new definition of distinguished trajectory that
generalises the concepts of fixed point and periodic orbit
to aperiodic dynamical
systems. This new definition is valid for identifying
distinguished trajectories with hyperbolic and non-hyperbolic types
of stability.
The definition is implemented numerically  and the procedure
consist in determining a path of {\it limit coordinates}.
It has been successfully applied to known examples of distinguished trajectories.
In the context of highly aperiodic realistic flows our definition
characterises {\it distinguished trajectories} in finite  time intervals,
and states that outside these intervals trajectories are no longer {\it distinguished}.

\end{abstract}

\baselineskip 20pt

\noindent Corresponding author: Ana M. Mancho,  Instituto de Ciencias Matem\'aticas,
Consejo Superior de Investigaciones Cient\'{\i}ficas,
Serrano 121, 28006, Madrid. e-mail: A.M.Mancho@imaff.cfmac.csic.es.
Phone: +34 91 5616800 ext 2408. Fax: +34 91
5854894.

\baselineskip 20pt \newpage

{\bf This paper  attempts to generalise the concepts of fixed point and periodic orbit
to time dependent aperiodic dynamical systems. Fixed points and periodic orbits are  keystones
for describing  solutions of autonomous
and time periodic dynamical systems, as the stable and unstable manifolds of these hyperbolic objects
 form  the basis of the
geometrical template organising the description of  the dynamical system.
The mathematical theory  of aperiodic dynamical systems is far from  complete. In this context, this work
deals with a general definition that encompasses the concepts of fixed point
and periodic orbit and which when  applied  to finite time and
aperiodic dynamical systems identifies special trajectories that play an
organising role in  the geometry of the flow. }

\section{Introduction}
 In recent years the theory of dynamical systems
has provided a useful framework for describing transport
in fluid  flows. Since the seminal work
by Aref \cite{Aref} on chaotic advection
much progress has been made both in theory and applications.
Dynamical systems techniques were first applied to Lagrangian transport in the context of two-dimensional, time-periodic flows
\cite{tpfl} and stationary 3D flows such
as the ABC flow \cite{3dsfl}.
More recently these techniques have been extended to describe aperiodic flows \cite{wiggins,malhotra,haller1} and finite
time-dependent flows, such
as those rising in geophysical applications \cite{physrep,nlpg}.
However, the mathematical theory for both aperiodic time-dependent flows and finite time aperiodic flows
 is far from being  completely developed.

For stationary flows the idea of {\it fixed point} is a key for describing
geometrically the solutions. Fixed points may be classified as hyperbolic or non-hyperbolic
depending on their stability  properties.
Stable and unstable manifolds
of hyperbolic fixed points organise the phase portraits of the flow away from the region close to the fixed points \cite{wigg,guck}.
These manifolds comprise respectively  the trajectories that
approach the fixed points as time tends to plus or minus infinite. As
they are formed of trajectories they act as barriers to transport
as particles cannot cross them without violating the
uniqueness of the solution. They are useful because
they allow qualitative predictions for the evolution of sets of initial conditions avoiding
explicit integration of
initial conditions on the whole domain.
Hyperbolic fixed points and their stable and unstable manifolds are the basic notions used for
the geometrical description of flows in autonomous dynamical system.

The concept of fixed point is extended to time periodic flows by means of the Poincar\'e map, as periodic orbits
with  period $T$ become fixed points of
the Poincar\'e map. For hyperbolic periodic orbits there exist also stable and unstable
manifolds that are geometric objects that organise the global dynamics.
Again they are respectively the sets of orbits asymptotically
approaching the periodic orbit as time tends to plus or minus infinity.

Aperiodic flows are still poorly understood, as theory that
is well established for autonomous or periodic flows does not apply  to them directly.
For instance  there exists  efforts in the mathematical community to extend the well known concept
of bifurcation for stationary flows  to non-autonomous systems \cite{langa1,langa2}.
To gain insight on the geometrical structure of  aperiodic flows, concepts such as Lyapunov exponents
are used, however these are  defined strictly on infinite time systems.
Realistic flows, like those arising  in geophysics or oceanography,
are not infinite time systems and for their description,  finite time versions of
the definition of Lyapunov exponents such as
 Finite Size Lyapunov Exponents (FSLE) \cite{aurell} and
Finite Time Lyapunov Exponents (FTLE)  \cite{haller,nese} are used.

Special trajectories, such as
detachment and reattachment points   \cite{haller2},
are observed  in  highly aperiodic  or turbulent flows.
In particular these separation trajectories occur on the boundaries
in simplified ocean models \cite{coulliete} and
also in technological applications in air foil design \cite{eis}.
Recent articles by Ide {\it et al} and Ju {\it et al} \cite{kayo,yu} referring
to these special trajectories
introduce the concept of {\it Distinguished Hyperbolic Trajectory } (DHT) which encompases
not only trajectories on the boundaries but also
special trajectories in the interior of the flow. DHT are
hyperbolic trajectories that, like hyperbolic  fixed points and periodic orbits,
have  stable and unstable manifolds that
are key for describing geometrically the solutions on the phase space.
This generalisation is an important step-forwards in the study of aperiodic flows, as it is a powerful tool
for describing transport in realistic oceanographic  flows \cite{physrep, nlpg, jpo,br}.
Distinguished hyperbolic trajectories as defined in \cite{kayo,yu}
are computed from hyperbolic instantaneous stagnation points (ISPs)
by means of an iterative procedure. If
instantaneous stagnation points bifurcate and do not persist for all times
the technique  developed in \cite{kayo,yu} cannot be applied in those time intervals, leaving many questions  unanswered,
such as what happens to the distinguished trajectories at those times,
for
distinguished hyperbolic trajectories
are trajectories, and as trajectories  exist at all times.
In fact  Ref. \cite{physrep} provides examples of  vector fields
with exact distinguished hyperbolic trajectories that exist on  time intervals without hyperbolic ISP.
Refs. \cite{physrep,nlpg}  discuss
 the impossibility of this technique for  tracking DHTs after  ISP
bifurcations and as a consequence the difficulties in establishing whether DHTs
obtained at different times are part of the same trajectory or not.

In this paper, following  ideas discussed in  \cite{physrep,kayo,yu},
we propose a new definition of  {\it Distinguished Trajectory } (DT)
which generalises the concepts of fixed point and periodic orbit to aperiodic flows.
We have taken the liberty of calling them {\it Distinguished} as in \cite{physrep,kayo,yu},
since although the definitions are not strictly
equivalent, it is found that  the studied hyperbolic  trajectories
are encompassed by both definitions.
We remark that our notion has the advantage over the method proposed in  \cite{kayo,yu}
that the DTs may be computed without the presence of hyperbolic instantaneous stagnation points.
Our definition  does not
depend on the  dimension $n$ of the space on which the vector field is defined and is valid both for hyperbolic and non-hyperbolic types of stabilities.
Non-hyperbolic DTs have not been studied in \cite{kayo,yu}, and in this sense our definition is broader
than that proposed there. In particular, we will show that exact non-hyperbolic periodic orbits
fall within the category of  distinguished trajectories.
Trajectories of this type could be of special interest for their applications in oceanography, as they are related to eddies and vortices.
Ocean eddies are well studied \cite{samelson}. Frequently they are
long lived, and water trapped inside can maintain its biogeochemical properties
for long time, being transported with the vortex. In steady horizontal velocity
fields, the presence of closed streamlines is the mathematical reason
for the isolation of the vortex core from the exterior fluid.
In two-dimensional, incompressible, time-periodic velocity fields the KAM tori
enclose {\it the core}, a region of bounded fluid particle motions   that do not mix with the
surrounding region \cite{wiggins}. But how to define an eddy from the Lagrangian point of
view in aperiodic flows? This is still an open question \cite{br,haller3} for which we  will discuss
new possibilities suggested by  the definitions given in this article.

The structure of the paper is as follows. Section 2 introduces the definition
of distinguished trajectory and explains its motivation in the context of 1D examples.
Section 3 explains the algorithm used to verify the applicability of our definition
of distinguished trajectories to the solutions of the periodically forced Duffing equation.
Details about technical issues arising from implementation of the definition are given.
Section 4 reports the results obtained in several  other 2D and 3D examples,
both periodic and non-periodic, hyperbolic and non-hyperbolic.
Section 5 discusses  results on realistic flows. Attention is paid to open questions on
distinguished trajectories  such as those mentioned above and pointed out in Refs. \cite{nlpg,physrep}.
Finally,  section 6 presents the conclusions.

\section{Distinguished trajectories: a definition}

We start by recalling the definition of distinguished hyperbolic trajectory
 provided  in \cite{kayo}.  Given the system:

\begin{eqnarray}
\frac{d {\bf x}}{dt} &=&{\bf D x} + g^{NL}({\bf x},t) \label{eq1} \, \, {\bf x} \in \, \mathbb{R}^n
\end{eqnarray}

Let ${\bf x}(t)$ be a trajectory of Eq. (\ref{eq1}) that
remains in a bounded region for all time. Then ${\bf x}(t)$ is said to
be a distinguished hyperbolic trajectory if:

  1. it is hyperbolic,

  2. there exists a neighbourhood $\mathcal{B}$ in the flow domain having
     the property that the DHT remains in $\mathcal{B}$ for all time, and
     all other trajectories starting in $\mathcal{B}$ leave $\mathcal{B}$ in finite time,
     as time evolves in either a positive or negative sense,

  3. it is not a hyperbolic trajectory contained in the chaotic
     invariant set created by the intersection of the stable and
     unstable manifolds of another hyperbolic trajectory.

\begin{remark}  If the data spans only a finite time interval, then the DHT
cannot be determined uniquely. Instead, there is a small
region in  ${\mathcal B}$ where the DHT can exist.
\end{remark}

In \cite{kayo} this setup is extended to general vector fields as follows. Coordinate
transformations are sought which put the system in the form of Eq. 1 and then
the previous definition is applied.

We give now our definition of  distinguished trajectory for a general vector field:

\begin{eqnarray}
\frac{d{\bf x}}{dt} &=& {\bf v}({\bf x},t), \label{adv1}  \, \,{\bf x} \in \, \mathbb{R}^n, \, t \in \, \mathbb{R}
\end{eqnarray}

We assume that ${\bf v}({\bf x},t)$ is  $C^r$ ($r \geq 1$) in ${\bf x}$ and continuous in $t$. This will allow for unique solutions
to exist, and also permit linearization, although linearization will not be used in our construction.

Before giving our definition of DT, we first need to introduce some notation and to make some definitions.
Let ${\bf x}(t)$ denote a trajectory of the system (\ref{adv1}) and
denote its components in $\mathbb{R}^n$
by $(x_1,x_2,...x_n)$.
For any initial condition ${\bf x^*}$   in an open set ${\mathcal B}\subset\mathbb{R}^n$,
consider the  function $M: {\mathcal B}\to\mathbb{R}$
\begin{equation}
M({\bf x^*})_{t^*,\tau}= \left(  \int^{t^*+\tau}_{t^*-\tau} \! \!\!\sqrt{\sum_{i=1}^n \left(\frac{d x_i(t)}{dt}\right)^2 }dt \right), \label{def:M}
\end{equation}
$M$ is the function that associates to each initial condition ${\bf x^*}$ in ${\mathcal B}$ the arclength of the trajectory that
passes through ${\bf x^*}$ at time $t^*$. The arclength of the trajectory is considered over its projection in the phase
space $(x_1,x_2,...x_n)$ and depends on $t^*$ and $\tau$.
As the function M is defined over an open set it does not necessarily attain a
minimum, but if it does, the minimum is denoted by  ${\rm min}(M({\bf x^*})_{t^*,\tau})$.

\begin{definition} ($\tau$-Distinguished trajectory).  A trajectory $\gamma(t)$
of Eq. (\ref{adv1}) is $\tau$-distinguished at time $t^*$  if there exists
an open set ${\mathcal B}$
around $\gamma(t^*)$ on which the defined function $M({\bf x^*})_{t^*,\tau}$ has a  minimum
  and
\begin{eqnarray}
{\rm min}(M({\bf x^*})_{t^*,\tau})=M( \gamma(t^*) )_{t^*,\tau}.
\end{eqnarray}\label{def:dt}
\end{definition}

\subsection{A discussion of the definition}
The
elements of the above definition deserve a detailed justification.
We illustrate our explanations with examples in 1D. First we consider an example taken from \cite{kayo,szeri}.
It is the linear one-dimensional non-autonomous dynamical system given by:
\begin{eqnarray}
\frac{dx}{dt}&=&-x+t. \label{eq1b}
\end{eqnarray}
For this example we consider the DHT reported  in \cite{kayo}, which is given by $x=t-1$. This is the particular
solution of the linear equation (\ref{eq1b}) towards which
all trajectories decay. The solution through the point $x^*$ at $t=0$ is  given by,
\begin{eqnarray}
x(t)=t-1+e^{-t} (x^*+1). \label{sol1}
\end{eqnarray}
Figure \ref{fig:1}a) displays several trajectories starting at
times ranging from $t = 0$ to $t = 4$ and figure \ref{fig:1}b) displays the same but starting at time
$t=-4$ (note that in this case part of the trajectories are out of the displayed domain). For each initial
condition the function $M$ provides the length of the projection of the trajectory
over the $x$-axis in the range of times $[-\tau, \tau]$.
Geometrically  it is clear that in this example  the function $M$ should have a minimum for a certain
$x$ value and that this value
 depends  on $\tau$. Ideally the minimum of $M$ should coincide
with the position of the DHT  at  $t=0$, however this would not be possible
if   in the definition of $M$
only positive times were considered, {\it i.e.} if the limits of the integration were $(0, \tau)$
the dashed trajectory in figure \ref{fig:1}a) would have a lower projection in positive times than the particular solution. An analogous problem would be encountered were only negative times considered,
that is if the limits of the integration
would have been $(-\tau,0)$. To determine precisely the position of the DHT at  $t=0$, both positive and negative times must be considered in the definition of $M$.

Figure 1b) confirms that with this choice the dashed trajectory cannot be distinguished, as it increases its projection
in negative times.
 Figure \ref{fig:2}a) displays   the
function $M(x^*)_{t=0,\tau}$  evaluated along the trajectories
(\ref{sol1}), for several $\tau$ values. Figure \ref{fig:2}b) displays the position of the minimum
of the  function $M_{t=0,\tau}$
as a function of $\tau$. These minima correspond to the positions of the $\tau$-distinguished
trajectories at $t=0$ and as $\tau$ increases
they  approximate the coordinate of the DHT at this time, which is at $x^*=-1$.
 The pair $(t^l, {\bf x}^l)$ formed by the time
at which $M$ is computed and the
value of the coordinate ${\bf x}^l$ to which
the minimum of the  function $M_{t^l,\tau}$ converges for increasing $\tau$
is called the {\it limit coordinates}.
The graphic \ref{fig:2}b)  illustrates the idea of approaching a point $(t_0, {\bf x}_0)$ of the distinguished trajectory by means
of the {\it limit coordinates}. In practice the convergence to the {\it limit coordinates} cannot be examined in the limit $\tau  \to \infty$,
either because it is impracticable in a numerical implementation,
or because in the large $\tau$ limit   errors accumulate, or simply
because the dynamical system is defined by a finite time data set. For these reasons
the convergence to the {\it limit coordinates} will be tested up to a finite $\tau$.

Figure \ref{fig:2}b)  raises
the question: what controls the
rate of the convergence of the minima of $M$ to the coordinates of the DHT?

It is hard to  answer  this question rigorously for a  vector field as general as in Eq. (\ref{adv1}).
However, some insight may be provided by particular examples.
 For instance  the system
\begin{eqnarray}
\frac{dx}{dt}&=&-2x+2t-1,
\end{eqnarray}
has the same DHT as (\ref{eq1b}). Its solution through the point $x^*$ at $t=0$ is  given by
\begin{eqnarray}
x(t)=t-1+e^{-2t} (x^*+1).\label{sol2}
\end{eqnarray}
Here the decay of the solution towards the DHT is faster due to the presence of the exponential
term $e^{-2t}$. Figure  \ref{fig:3} shows that in this case  the
rate of the convergence of  the minima of $M$ towards the coordinates of the  DHT
at time $t=0$ is also faster than before. However
there exist systems in which  the exponential decay of the solution is not a determining
 factor affecting the  rate of the  convergence of the minima of $M$
to the coordinates of the  DHT. For instance, in autonomous systems fixed points are the DTs, and clearly
they are minimizers of $M$ for any $\tau>0$ whatever is  the
exponential rate of growth or decay    of the nearby solution.

In these examples the function $M$ has a unique minimum, but as we will see the situation
will  not always be
so simple when nonlinearities are involved in the vector field. Also it is important to notice that
the function $M$ obtained at different $\tau$ values has been  used to obtain the
{\it limit coordinates} $(t^l_{0},{\bf x}^l_0)$ and that these approach the ${\bf x}_{0}$ coordinate of DHT {\it at a given
time $t_0$} (here $t_0=t_0^l)$. Once this is obtained,  approaching the DHT at later times $t_k=t_0+ k \Delta t$ would require
applying the same procedure to get the {\it limit coordinates} $(t^l_{k},{\bf x}^l_k)$.
We remark here that the proposed algorithm does not ensure that the set of  {\it limit coordinates}  $(t^l_k, {\bf x}^l_k)$
are in fact part of a trajectory. Later we will see that in practice, in many examples these
points approach a true trajectory,  however in realistic aperiodic flows this has to be
verified {\it a posteriori}. These considerations lead us to the definition of a {\it Distinguished Trajectory}.
\begin{definition} (Distinguished trajectory). A trajectory $\gamma(t)$ is said to be Distinguished with accuracy $\epsilon$ ($\,0\leq \epsilon$ )
in a time interval $[t_0, t_N]$ if there exists  a continuous path  of {\it limit coordinates} $(t^l,{\bf x}^l)$  where
$t^l \in [t_0, t_N]$, such that,
\begin{eqnarray}
||\gamma(t^l)-{\bf x}^l(t^l)|| \leq \epsilon, \,\, \forall t^l \in [t_0, t_N].
\end{eqnarray}\label{def:dtdef}
Here $|| \cdot ||$ represents the distance defined by
\[||{\bf a}-{\bf b} ||=\sqrt{\sum_{i=1}^{n} (a_i-b_i)^2} \,\,\,\, \,{ with} \,\,\,\,{\bf a}, {\bf b} \in \mathbb{R}^n.\]
\end{definition}
In the numerical exploration of this definition
we will replace the continuous path  of {\it limit coordinates} $(t^l,{\bf x}^l)$ and the continuous trajectory $\gamma(t)$ by
discrete representations  $(t^l_{k},{\bf x}^l_k)$ and $\gamma(t^l_k)$ where $t_0\leq t^l_k \leq t_N$.
By definition \ref{def:dtdef}  any trajectory is
distinguished for sufficiently large  $\epsilon$, however the interesting distinguished trajectories
 are those for which $\epsilon$ is close to zero, which means
it is of the order of the accuracy in which $\gamma(t^l_k)$ and ${\bf x}^l(t^l_k)$ are numerically determined,
or zero, if an exact expression is known for both.

Underlying definitions 2 and 3 is the geometrical idea that distinguished trajectories,
which act as organising centres of the flow in phase space, are those
that "move less" (in a certain sense) than other nearby trajectories.
 This property of ``moving less'' is
satisfied by minima of the  function $M$ as it measures
the length of the displacement in phase space of a trajectory forwards and backwards in time.
In fact this property is related somehow to property (2) of the definition
provided in \cite{kayo} and presented at the beginning of Section 2, as the trajectory that ``moves least''
is not  expected to leave  the neighbourhood $\mathcal{B}$.

Definitions \ref{def:dt} and \ref{def:dtdef} are made for a general dynamical system in any dimension $n$.
The purpose of this paper is the exploration of these definitions, but more in an illustrative than demonstrative way,
as it is impossible to provide examples
for every possible $n$, and  one cannot deal with every possible example at a given $n$.
Even if one wants to provide a rigorous formal proof that  the definition  recovers specific trajectories
such as periodic orbits (it is not obvious that in general they have to satisfy our definition),
this has to be done with some further hypotheses
on the vector field and  proofs
will not be valid beyond the assumed hypotheses. Therefore we restrict the
 discussion to dimensions up to 3, as these are the  dimensions important for
geophysical flows, which are what originally
motivated the definition. However  it is sensible to make  the same definition
for any  dimension $n$, as it is clear that it works
for autonomous systems of any dimension.
Fixed points are  the kind of trajectory  expected to be recovered
by the definition and they do not move at all in the phase space.
For these  $M=0$, while $M>0$ for any other trajectory in the neighbourhood which is not a fixed point.

We conclude this section with some remarks. First, it is not guaranteed {\it a priori} that for an arbitrary
vector field, satisfying only some  rather general hypotheses such as those of Eq. (\ref{adv1}),
the function $M$ will have  a minimum, however this is not a problem from the point of view
of  the definition.
For instance the same thing happens for general non-linear autonomous systems.
 In these systems fixed points are perfectly
defined although one does not know {\it a priori} if such points  exist for arbitrary examples.
If they exist, it is possible to find them by either solving the nonlinear equation ${\bf v}({\bf x})=0$,
or by applying  definitions
\ref{def:dt} and \ref{def:dtdef}.
In the same way one does not know  {\it a priori} if distinguished trajectories  exist for
a general vector field
however if they exist they can been found with the tools proposed in this article.
Second, even if a path of {\it limit coordinates} is found it
is not guaranteed that it will be a trajectory, although if that is the case then
from definition \ref{def:dtdef} follows  that this trajectory is distinguished.
Third, one might think that if {\it limit coordinates} are found
at $t_0$  that  approach with great accuracy a point of
{\it an existing DT}, then the iterative procedure described above for finding
a set  of {\it limit coordinates} $(t^l_{k},{\bf x}^l_k)$
approaching the DT at later times is an unnecessary computational effort,
as those coordinates could have been equally well obtained
by integrating forwards the initial data. However there exist examples such as  a hyperbolic
DT in dimension greater than one with at least 1D unstable manifold, that cannot be integrated like this, as
 the integrated trajectory will
eventually leave the neighbourhood of the DT through the unstable manifold no matter how small the initial error is.
In summary  the proposed methodology
based on {\it limit coordinates} provides a systematic way of finding DT, which
can be elusive and difficult to obtain.
We will discuss these issues in detail in later sections.

\section{A numerical algorithm}

In this section we propose an algorithm for computing
a path of {\it limit coordinates} in a time interval, and
we verify that it is close to a DT of a known example.

For this purpose  we calculate, at increasing $\tau$ values, the minimum of the function $M_{t=0,\tau}({\bf x})$ for ${\bf x}$ in  an open set in $\mathbb{R}^n$.
The method is  illustrated in a  2D case, the periodically forced Duffing equation
\begin{eqnarray}
\dot{x}&=&y,\nonumber \\
\dot{y}&=&x-x^3+\varepsilon\sin (t),\label{eq:duffing}
\end{eqnarray}
where $\varepsilon$ is a small parameter. The hyperbolic fixed point of the unperturbed  autonomous system (i.e., $\varepsilon=0$)
is at the origin ${\bf x}= (0,0)$. For small $\varepsilon$, it is
possible to compute by perturbation theory (see \cite{physd}),
the following periodic trajectory which stays close to the origin:

\begin{equation}\label{eq:dhtduffing}
{\bf x}_{DHT}(t)=-\frac{\varepsilon}{2} \binom{\sin t}{\cos t}-\frac{\varepsilon^3}{40}
\binom{2\sin^3 t+\frac{3}{2}\sin t\cos^2 t}{\frac{3}{2}\cos^3 t+3\sin^2t\cos t}
+\mathcal{O}(\varepsilon^5).
\end{equation}

For $\varepsilon=0.1$,  Eq. (\ref{eq:dhtduffing})
is accurate   up to the fifth digit. This trajectory  is  identified as {\it Distinguished} in Ref. \cite{kayo}, for this reason we
have labelled it a $DHT$.
Substituting the expression,
\begin{equation}
{\bf x}=(x,y)={\bf x}_{DHT}(t)+\left(\xi_1,\xi_2 \right)
\end{equation}
into Eq. (\ref{eq:duffing}) and by dropping the nonlinear terms one finds that the linearized
equations have two linearly independent solutions in terms of  which the time
evolution of the components ${(\xi_1, \xi_2)}$ is:
\begin{equation}
\left(\xi_1,\xi_2 \right) = \alpha \,{e}^{t} \left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2} \end{array}\right)+\beta \, {e}^{-t} \left(\begin{array}{r}-1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)+\mathcal{O}(\varepsilon^2).\label{timeev}
\end{equation}
Eq. (\ref{timeev}) confirms the hyperbolicity of the solution (\ref{eq:dhtduffing}).

This explicit expression for the distinguished hyperbolic trajectory is
a benchmark for testing the utility of our definition.
The procedure starts by determining the coordinates of ${\bf x}_{DHT}$ at time $t=0$.
We consider the open set ${\mathcal D}\subset \mathbb{R}^2$, defined by
 ${\mathcal D}=(-0.2,0.2)\times(-0.2,0.2)$ and in the function $M_{t=0,\tau}({\bf x})$ we take $\tau$
to be $2$.
Figure \ref{fig:duffinggrid} displays a contour plot of $M_{t=0,\tau=2}({\bf x})$ which has a minimum  at ${\bf x}=(0, -5.7057  \cdot 10^{-2})$. $M_{t=0,\tau=2}({\bf x})$ quantifies displacements of particles in  phase space, and its minimum corresponds to
the initial condition that
``moves less'' over the $\tau$ interval $[-2, 2]$.
As noted in the previous section,  when the
 value of $\tau$ is increased,
the position of the minimum gets closer and closer to  the coordinates of the DHT at $t=0$.
 Figure \ref{fig:duffingtau} shows contour plots of the function $M$ for several $\tau$ values.
Fig.  \ref{fig:duffingtau}a) displays a typical hyperbolic structure for $M$
 for $\tau=5$ where the directions of the stable and unstable manifolds are easily recognised.
In Fig.  \ref{fig:duffingtau}a) the function $M$ has a unique minimum  at  ${\bf x}=(0, -4.979 \cdot 10^{-2})$
while in Fig.  \ref{fig:duffingtau}b) there  appear several minima for $\tau=10$.
The global minimum in this picture corresponds to ${\bf x}=(0, -5.0042565261 \cdot 10^{-2})$.
Figure \ref{fig:duffingtray}a)  compares the $x$-coordinate of ${\bf x}_{DHT}$
as a function of time with  trajectories having initial conditions
at the global minima of  $M_{t=0, \tau=2}$ and $M_{t=0, \tau=10}$.
Taking  as initial condition the global minimum of $M_{t=0, \tau}$  for $\tau=10$
provides a trajectory that  stays  close to
 ${\bf x}_{DHT}$ for a longer time interval than  for $\tau=2$, which
confirms  that larger $\tau$-values more closely approach the coordinates of the DHT.
Fig.  \ref{fig:duffingtau}c) displays the contour plot of $M_{t=0,\tau=50}({\bf x})$. Its
global minimum is  at ${\bf x}=(0,-5.0037606418 \cdot 10^{-2})$.
The associated trajectory depicted in Fig. \ref{fig:duffingtray}b) shows that
this initial condition tracks the DHT for a longer time interval than those obtained for $\tau=2$ and 10,
however the figure shows that  the integration of the DHT
in  $(-50, 50)$ is not possible. In fact the associated
trajectory stays close to the DHT  only in  the time interval $(-20, 20)$. This confirms that results
obtained for $\tau=50$ are the same as those obtained for $\tau=20$.
In practice for a finite precision  numerical scheme,
such as
a 5th order Runge Kutta used here, the approach to the DHT has
an upper bound depending on $\tau$. This occurs because the stable and unstable manifolds of the  hyperbolic
trajectory magnify any initial error in either negative  or positive  time and
beyond this $\tau$-limit numerical  errors dominate. The convergence towards
the DHT is confirmed in Fig. \ref{fig:duffingxy} which
displays the evolution of the coordinates $x$ and $y$ of the global minimum of $M$ as a function of the parameter $\tau$.

New minima appearing in Figs.  \ref{fig:duffingtau}b) and c)
relate to the existence of different $\tau$-distinguished trajectories.
As illustrated in Fig. \ref{fig:duffingtray}b), they correspond to trajectories which stay close to
${\bf x}_{DHT}$ in a small time range contained in the interval $-\tau < t < \tau$,
but which later fly apart from the DHT.

We now describe a numerical scheme to compute a path of limit coordinates.
The algorithm has the following steps:

\begin{enumerate}
\item Step 1. Discretize the domain  ${\mathcal D}$ at the initial time $t=t_0$ at which
one wishes to compute a DT. For instance, the grid size of this domain in Fig.  \ref{fig:duffinggrid}
is $101 \times 101 $. The function $M$ is evaluated  at each grid point for a given  $\tau_0$.
\item Step 2. Search for the local minima of  $M_{t_0,\tau_0}$ in the interior of the
grid. These minima approach   the coordinates of $\tau_0$-distinguished trajectories within the accuracy of the grid.
In what follows we restrict our description to the case of a unique minimum,
as this simplifies the description;
the procedure is easily generalised to the case of multiple minima.
\item  Step 3. Improve the approach of the coordinates of the $\tau_0$-distinguished trajectory up to precision $\delta$.
For this purpose build up a
 $3^n$ grid
centred on the candidate point provided by step 2, (for the 2D case this is a $3\times 3$ grid as Fig. \ref{fig:epsilongrid} illustrates),
setting the  distance between nodes equal to $\delta$.
Then evaluate $M_{t_0,\tau_0}$ at the points of the $\delta$-grid. If the minimum  of $M_{t_0,\tau_0}$ is in the interior
of the grid, then the coordinates of the $\tau_0$-distinguished trajectory are known to within $\delta$ accuracy. Otherwise
the $\delta$-grid must be rebuilt centred on the boundary point where the minimum has been located, and
$M_{t_0,\tau_0}$ must be re-evaluated in the new $\delta$-grid. This procedure stops when the minimum
of $M_{t_0,\tau_0}$ is in the interior of the grid.
\item  Step 4. Computing the {\it limit coordinates} at time $t_0$.
Define a sequence of increasing $\tau$-values as follows: $\tau_1=\tau_0+ \Delta \tau$ and $\tau_2=\tau_0+ 2 \Delta \tau$.
Then evaluate $M_{t_0,\tau_0}$, $M_{t_0,\tau_1}$ and $M_{t_0,\tau_2}$
on the $\delta$-grid. If  the minimum is at  an interior position for the three cases,
then we consider that {\it limit coordinates} have been found within  $\delta$  accuracy.
We note that this is a necessary but not sufficient condition as one does not know {\it a priori}
the convergence  rate to the distinguished trajectory. Although this criterion could be strengthened,
it has been tested and found to be adequated for the examples explained in subsequent sections.
If the condition  defined above of having a minimum  at  an interior position for the sequence of $\tau$-values
is not satisfied then, after  replacing $\tau_0$ by $\tau_1$, we return to step 3 and then to step 4.
The loop between steps 3 and 4 is stopped when the condition of step 4 is satisfied for some $\tau_k$.
\item Step 5. Compute the {\it limit coordinates} at time $t_1=t_0 + \Delta t$. Once the {\it limit coordinates}
have been approached
at time $t_0$, they are  integrated forward  numerically up to time $t_1$. If the {\it limit coordinates} converge to a
hyperbolic  DT
with an unstable manifold, the position ${\bf x}(t_1)$ obtained should  deviate from the position
of the DT at time $t_1$. In order to correct this, the procedure described above is repeated from step 3 onwards.
For that purpose in the definition of $M$, $t_0$ is replaced by $t_1$ and the $\tau$-value
is reset to $\tau_0$. The coordinates ${\bf x}(t_1)$ are the first
 approximation to the  $\tau_0$-distinguished trajectory at time $t_1$. Once the {\it limit coordinates}
are found
for time $t_1$ it is possible to repeat the procedure  to locate them at successive times
 $t_2,t_3,\ldots,t_N$.
\end{enumerate}

The  algorithm  requires as inputs: an explicit expression for
the dynamical system (\ref{adv1}); the definition of the domain $\mathcal{D} \subset \mathbb{R}^n$; the initial and final
times $t_0$, $t_N$ at which DTs are required,
and the time step $\Delta t$ for intermediate times; the initial $\tau_0$
and the increment $\Delta \tau$; the precision $\delta$. As an output the algorithm gives
a path of limit coordinates  at the selected times $t_k$.

Next we discuss in more detail some technical aspects related to the implementation
of the above algorithm. Steps 1 and 3 require evaluating
$M_{t_0,\tau_0}$ as defined in Eq. (\ref{def:M}). We explain how this is done for the contour plots displayed
in Figs.  \ref{fig:duffinggrid} and \ref{fig:duffingtau}, which  refer
to the system (\ref{eq:duffing}) at $t_0=0$.

Fig. \ref{tra1} shows a schematic projection onto the $\mathbb{R}^2$ plane of a possible trajectory
${\bf x}(t)$ of the system  from $-\tau$ to $\tau$.
As it was obtained  numerically, only a finite number of points ($L$) appear. This picture suggests the following  discrete version of Eq. (\ref{def:M}) for $M$:
\begin{equation}
M({\bf x})_{0,\tau}= \sum_{j=1}^{L-1} \left(  \int^{p_f}_{p_i} \! \!\!\sqrt{\left(\frac{d x_j(p)}{dp}\right)^2+\left(\frac{d y_j(p)}{dp}\right)^2 }dp \right),
\end{equation}
where the functions $x_j(p)$ and $y_j(p)$ represent a curve interpolation parametrized by $p$, and the integral
\begin{equation}
\int^{p_f}_{p_i} \! \!\!\sqrt{\left(\frac{d x_j(p)}{dp}\right)^2+\left(\frac{d y_j(p)}{dp}\right)^2 }dp
\end{equation}
is computed  numerically. In our case we use the Romberges method (see \cite{nr}) of  order $2K$ with $K=5$.
It is clear that the accuracy of the evaluation of $M$ will depend on the number of points
 on the trajectory  $L$, which is controlled by the size of the time step, $h$, of the integrator (a 5th order Runge Kutta method)
and on the interpolation scheme between points.
Two interpolation methods are compared in tables (\ref{tab:erroresl}) and (\ref{tab:erroresd}).
Results in table (\ref{tab:erroresl}) are obtained with linear interpolation between nodes. Results in
table (\ref{tab:erroresd}) correspond to the interpolation method used by Dritschel \cite{dr} in the context of contour dynamics,
 which has been successfully applied in \cite{physd} to the computation of invariant manifolds
for aperiodic flows.
This method interpolates a piece of the curve in Fig. \ref{tra1} between consecutive nodes as follows:
\begin{equation}
{\bf x}_j (p) = {\bf x}_j + p {\bf t}_j + \eta_j (p){\bf n}_j \label{eq:intdri}
\end{equation}
for $p_i=0 \leq  p \leq  p_f=1$ with ${\bf x}_j (0) = {\bf x}_j$ and ${\bf x}_j (1) = {\bf x}_{j+1}$, where:
\begin{eqnarray}
{\bf t}_j &=& (a_j , b_j ) = {\bf x}_{j+1} - {\bf x}_{j},  \,\,\, {\bf t}_j  \in \mathbb{R}^2\\
{\bf n}_j &=& (-b_j , a_j ),   \,\,\,  {\bf n}_j \in  \mathbb{R}^2\\
\eta_j (p) &=& \mu_j p + \beta_j p^2 + \gamma_j p^3, \,\,\,  \eta_j \in    \mathbb{R}.
\end{eqnarray}
The cubic interpolation coefficients $\mu_j$, $\beta_j$ and $\gamma_j$ are
\[
 \mu_j = - \frac{1}{3} d_j \kappa_j - \frac{1}{6} d_j \kappa_{j+1}, \,\,\, \beta_j = 2 d_j \kappa_j, \,\,\, \gamma_j = \frac{1}{6}d_j
(\kappa_{j+1} - \kappa_j ),
\]
where $d_j = |{\bf x}_{j+1} - {\bf x}_j |$ and
\[
\kappa_j = 2 \frac{a_{j-1} b_j - b_{j-1} a_j}{|d^2_{j-1} {\bf t}_j + d^2_j {\bf t}_j-1 |}                                         \]
is the local curvature defined by the circle through the three points, $x_{j-1}$, $x_{j}$, and $x_{j+1}$.

Tables (\ref{tab:erroresl}) and (\ref{tab:erroresd}) show the errors in
the computed lengths of the ellipses for different  ratios of
major to minor axis.  The reference length is that obtained with GNU Octave
version 2.1.73, as it provides 16 correct digits for the known circumference.
The tables confirm that the Dritschel's method is superior to linear interpolation and
it is the one used to
compute the function $M$. In the trajectory from $-\tau$ to $\tau$ the number of points $L$ is determined by the time step size
of the Runge Kutta method which is set to $10^{-2}$.

Another important element  of the algorithm needing
 discussion is the value of the input parameters, in particular of $\tau_0$ and
$\Delta \tau$. It is clear from Fig. \ref{fig:duffingtau} that large $\tau$ values are not convenient as they
increase the roughness of the function
$M$ and  several local minima may appear in the neighbourhood of a DHT  that correspond
to trajectories that stay close to it
for some time. On the other hand it is clear that sufficiently large  $\tau$ values are required to fix the coordinates of the DHT to within prescribed accuracy. Combining these observations suggests the
use of relatively small values for
the initial $\tau_0$.
In the example above $\tau_0=2$, provides, as a starting point, a smooth $M$  as that of Fig. \ref{fig:duffinggrid}.
The increments should not be large. In practice we have chosen $\Delta \tau=\tau_0/2$. This prevents
from stepping to a  too rough $M$ before getting close enough to the sought after DHT. Some of the  local
minima appearing in Fig. \ref{fig:duffingtau}b) are just apparent and disappear with a more refined grid.
However, as already observed,
others  belong to  true $\tau$-distinguished trajectories, which are secondary and can be  avoided
if  the increment of the $\tau$-values is conveniently small. These choices are found to be appropriate
for determining with great accuracy the DHT in (\ref{eq:dhtduffing}) by means of a path of {\it limit coordinates}.
Figure \ref{fig:dhtduffing}a) represents both the
analytical DHT and the numerical limit coordinates
and Figure \ref{fig:dhtduffing}b) displays the distance  between the
exact and the numerical approach, confirming that the DHT in (\ref{eq:dhtduffing}) is also a DT
in the sense of our definition
\ref{def:dtdef} with
accuracy $\epsilon=10^{-6}$. Other parameters in the
algorithm are: $\delta=10^{-6}$, step size in the Runge Kutta method $ h=10^{-2}$,  $t_0=0$, $t_N= 6$, and $\Delta t=0.01$.
To locate the DHT with accuracy $\delta=10^{-6}$ requires increasing values of $\tau$  up to  15, which is
near  the limit of the integration method. Figure \ref{fig:dufftau}
shows the maximum $\tau$ required at each $t_k$.

\section{Applications to exact examples}

In this section we apply the algorithm explained in the previous section to selected examples.

\subsection{ A non-hyperbolic distinguished trajectory}

The unperturbed autonomous system (\ref{eq:duffing}) obtained with $\varepsilon=0$ has  non-hyperbolic
fixed points at $(-1,0)$ and $(1,0)$. Obviously these fixed points correspond to DTs
which are also $\tau$-distinguished trajectories for all $\tau>0$.
For the periodically forced system (\ref{eq:duffing}) with
small $\varepsilon$ it is possible  using perturbation theory
to find periodic solutions close to these fixed points in a manner similar to the analysis
of the hyperbolic example  made in the previous section.
For instance close to the point $(1,0)$ we find the periodic trajectory:
\begin{equation}\label{eq:detduffing}
x_{DET}(t)=-\binom{1}{0}+\varepsilon \binom{\sin t}{\cos t}+3\varepsilon^2
\binom{\frac{1}{2}\cos^2 t}{-\sin t\cos t}
+\mathcal{O}(\varepsilon^3).
\end{equation}
This solution has not been considered {\it Distinguished} in previous works \cite{kayo,yu}, as these have been focused
on hyperbolic trajectories and this solution, as is proved next,  is not hyperbolic.
However, in anticipation of its having the {\it Distinguished}
property, we have labelled it $DET$ for two reasons.
 One is that it is periodic, and we expect
periodic orbits to be distinguished, and second is that it is in clear correspondence to the elliptic fixed point (-1,0)
in the case $\varepsilon=0$, and fixed points are DTs.

To determine the  stability of  (\ref{eq:detduffing}) we proceed as before, by substituting the expression
\begin{equation}
{\bf x}=(x,y)={\bf x}_{DET}(t)+\left(\xi_1,\xi_2 \right)
\end{equation}
into Eq. (\ref{eq:duffing}). We find that the  linearized system at order $\varepsilon^0$ is:
\begin{eqnarray}
\frac{d \xi_1}{d t}&=&\xi_2\\
\frac{d \xi_2}{d t}&=&-2\xi_1.
\end{eqnarray}
Therefore the linearized flow around ${\bf x}_{DET}(t)$ evolves according to:
\begin{equation}
\left(\xi_1,\xi_2 \right) = \alpha \,{e}^{i\sqrt{2}t} \left(\begin{array}{c}1/\sqrt{3}\\ i\sqrt{2/3} \end{array}\right)+\alpha^* \, {e}^{-i\sqrt{2}t} \left(\begin{array}{r}1/\sqrt{3}\\ -i\sqrt{2/3}\end{array}\right)+\mathcal{O}(\varepsilon),
\end{equation}
which clearly is not hyperbolic. Here $\alpha$ and $\alpha^*$ are complex conjugate numbers.

We apply our algorithm to determine the {\it limit coordinates}
approaching (\ref{eq:detduffing}), as we want to verify whether
 definition~\ref{def:dtdef}
also  works for time-dependent non-hyperbolic solutions. The following input
 is considered: ${\mathcal D}=(-1.2,-0.8)\times(-0.2,0.2)$, $\tau_0=2$, $\Delta \tau=1$, $\delta=10^{-4}$,
$t_0=0$, $t_N=6$, and time step $10^{-2}$ for the Runge-Kutta integrator.
We note that the accuracy $\delta$ is not
as demanding as before, since now the exact ${\bf x}_{DET}$ for $\varepsilon=0.1$
is only accurate up to the third digit.
Figure \ref{fig:duffingetau} shows a rather different structure for the function $M$. An important
feature is the smoothness of $M$ close to the DET even for large $\tau$.
In figure \ref{fig:duffingetau}b)  the differences
between the  rather flat region around the position of the DT given by (\ref{eq:detduffing}), which appears in the dark tone,
and the roughness of the outer part are remarkable. The irregularity of this region suggests that inside it nearby trajectories follow
rather different paths as  happens for chaotic motions, while the regularity of the central core suggests the
existence of trapped trajectories circling around  the  DET. From
this perspective the function $M$ for large $\tau$
seems a useful tool for fixing
the boundaries of a Lagrangian eddy, different to the methods proposed in \cite{br,haller3}.

Figure \ref{fig:duffingexy} shows the rate of convergence to the global minimum of $M$  in the domain  ${\mathcal D}$
as a function of $\tau$. The convergence towards the coordinates of the DT  is oscillatory and rather
slow since $\tau$ values up to $600$ are required. A slight difference between the exact coordinates of the DT and
the numerically computed limit coordinates is evident, however we note that these differences are
consistent with  the precision to which the exact DT is known, which is only  to the third digit.
Figure  \ref{fig:detduffing}, and more specifically  figure  \ref{fig:detduffing}b),
confirms that the exact expression in Eq. (\ref{eq:detduffing})
is in fact a distinguished trajectory according to our definition~\ref{def:dtdef} with
accuracy $\epsilon=4 \cdot 10^{-3}$.

Figure \ref{fig:duffingetray}a) shows a forward and backward integration along the time interval $(-50, 50)$
taking as  initial data the limit coordinates supplied by our algorithm
at time $t_0=0$, and compares it with the exact solution of the DT.
From figure \ref{fig:duffingetray}b) it can be seen that this trajectory evolves  close to the exact
solution in the entire time range.
This result shows that contrary to what happens near hyperbolic trajectories,
near non hyperbolic trajectories,  small error does not amplify and as consequence,
once a DT {\it is known to exist} it could have been computed simply by integrating forwards and backwards
the limit coordinates found at a given time $t_k$. However one needs to be careful here, as a trajectory
is not necessarily distinguished at all times, and for it to be properly called
distinguished, it  should be verified that it stays close to the
{\it limit coordinates} in the whole time interval, and therefore one cannot avoid computing {\it limit coordinates}
along the time interval in this case either. We will return to this point in the next section.

\subsection{ The rotating Duffing equation}

Next we analyse the aperiodic hyperbolic distinguished
trajectory of a system already studied in \cite{physd},  the rotating Duffing equation:
\begin{eqnarray}
\left(\begin{array}{c}\dot{\eta_1}\\ \dot{\eta_2} \end{array}\right) &=& \left(\begin{array}{cc}\sin 2 \omega t & \cos 2\omega t +\omega\\
\cos 2\omega t -\omega & -\sin 2 \omega t
 \end{array}\right)\left(\begin{array}{c}\eta_1\\ \eta_2 \end{array}\right) \nonumber\\
&+&( \varepsilon \sin t - [ \cos \omega t
\eta_1 - \sin \omega t \eta_2 ]^3 )\left(\begin{array}{c}\sin \omega t \\ \cos \omega t \end{array}\right).
\end{eqnarray}
This Duffing equation is quasi-periodic in time when the rotation rate $\omega$ is irrational. It is obtained from
the system (\ref{eq:duffing}) by applying the rotation  ${\bf x} = R(t) \boldsymbol{\eta} $, where
\begin{eqnarray}
R(t)=\left(\begin{array}{cc}\cos \omega t & -\sin \omega t \\
\sin \omega t  & \cos \omega t
 \end{array}\right).
\end{eqnarray}
The DHT can also be obtained through the coordinate transformation:
\begin{equation}
{\boldsymbol{\eta}}_{DHT} (t) = R(t)^{-1} {\bf x}_{DHT} (t). \label{dhtrduff}
\end{equation}
Figure \ref{fig:dhtrduffing}, and in particular Fig. \ref{fig:dhtrduffing}b), confirms that the DHT (\ref{dhtrduff})
is also a DHT according to our definition ~\ref{def:dtdef} with accuracy $\epsilon=4 \cdot 10^{-6}$.

\subsection{ A 3D extension of the Duffing equation}

In this section we apply our definitions to an example in higher dimension. In particular we
consider a 3D extension of the Duffing equation:
\begin{eqnarray}
\dot{x}&=&y,\nonumber \\
\dot{y}&=&x-x^3+\varepsilon\sin (t),\label{eq:duffing3d}\\
\dot{z}&=&z+\varepsilon\sin (t).\nonumber
\end{eqnarray}
The hyperbolic fixed point of the unperturbed  autonomous system (i.e., $\varepsilon=0$)
is at the origin ${\bf x}= (0,0,0)$. The solution   for small $\varepsilon$ becomes:
\begin{eqnarray}

{\bf x}_{DHT}(t)&=&-\frac{\varepsilon}{2} \left(\begin{array}{c} \sin t\\ \cos t \\ {\small \cos t - \sin t}\end{array}\right)\nonumber\\&-&\frac{\varepsilon^3}{40}
 \left(\begin{array}{c} {2\sin^3 t+\frac{3}{2}\sin t\cos^2 t}\\{\frac{3}{2}\cos^3 t+3\sin^2t\cos t}\\ 0\end{array}\right)\label{eq:dhtduffing3d}
+\mathcal{O}(\varepsilon^5).

\end{eqnarray}

The numerical scheme explained in Section 3 is easily adapted to higher dimensions.
However some changes must be made.
The computation of $M$ requires approximating lengths of  trajectories which in 3D
needs an  interpolation scheme different to that
of Eq. (\ref{eq:intdri}), which is only valid in $\mathbb{R}^2$. We consider the linear interpolation instead.
This interpolation evaluates the function $M$ satisfactorily if trajectories are represented by  a large number of points.
This is achieved by using a Runge-Kutta  method with time step  $h=10^{-4}$.
Figure \ref{fig:duffingtau3d} indicates the evolution of  coordinates associated with the minimum of $M$
as a function of $\tau$ (solid line). The dashed line corresponds to the exact perturbative solution. There is evident
a clear  convergence towards the exact position although there is a significant jump in the asymptotic behaviour
beyond $\tau \sim 50$. This jump  is due to round off errors in the determination of $M$ for large $\tau$.
The third equation in (\ref{eq:duffing3d}) is just a linear equation and for this reason
 solutions which are in the neighbourhood of the DHT have $z$-coordinate
 growing exponentially  in backwards time. Thus for large $\tau$ values, the evaluation of  $M$ is made along  very long trajectories in the $z$-coordinate,
which are underrepresented by points   sampled every $h=10^{-4}$ (see table I) and where lengths
are badly calculated  by adding up  very small and very large (and inaccurate) numbers.
In spite of  this, figure \ref{fig:duffing3d} confirms that  the exact distinguished
trajectory can be accurately obtained with our methodology and that for $\tau<50$ errors are within the expected margin.
The remaining input parameters used  in Figs. \ref{fig:duffingtau3d} and \ref{fig:duffing3d}
are: ${\mathcal D}=(-0.2,0.2)\times(-0.2,0.2)\times(-0.2,0.2)$, $\tau_0=2$, $\Delta \tau=1$, $\delta=10^{-6}$,
step size $ h=10^{-4}$ in the Runge Kutta method, $t_0=0$, $t_N= 6$, and $ \Delta t=10^{-2}$.

As we explain next, the computational demands made by this example are considerably larger than
they were for the previously considered 2D examples.
When  determining a DT, most of the CPU time is spent computing the value of $M$ on
 the $\delta$-grid  displayed in Fig. \ref{fig:epsilongrid}. The number of neighbours
of the interior point grows with the dimension $n$ as $3^n$, therefore
when  the problem increases its dimension from $n$ to $n+1$, the computational demands are multiplied by 3.
Another factor that contributes to increased computational time is the decrease
of the Runge Kutta time step $ h$
in the evaluation of trajectories on the $\delta$-grid. This increases the number of points in the trajectory
(and therefore the number of  operations) with respect to the previous
Dristchel approach by a factor $100$. This factor  is partially balanced by the fact that for the same number of points
the arclength is computed more rapidly with the linear than with the Dristchel interpolation.

\section{Application to vector fields defined as finite time data sets}

In this section we explore definition \ref{def:dtdef} for a highly aperiodic 2D  flow in which
the vector field is defined as a finite time data set. In particular we consider the output of
a  quasigeostrophic wind-driven double gyre  model in a regime already studied in \cite{physrep,nlpg}.
Details of this model may be found in  \cite{nlpg, coulliete}. Fig.~\ref{fig:sfqgn} shows a typical
output for the streamfunction provided by this model.
The velocity data set is obtained on a 1000 km $\times$ 2000 km rectangular domain and spans 4000 days. This interval is considered for a fluid started from rest and allowed to spin for 25000 days.
Free slip  conditions are considered for the velocities on the boundaries
 and the wind stress curl is 0.32 dyn/cm$^2$.
The equations of motion for this system are given by:
\begin{eqnarray}
\dot{x}&=&v_x(x,y,t)=-\frac{\partial \psi}{\partial y},\\
\dot{y}&=&v_y(x,y,t)=\frac{\partial \psi}{\partial x},
\end{eqnarray}
and the variables $x$ and $y$ are in the rescaled domain $[0, 1]\times [0, 2]$. Here the  velocity fields $v_x$ and $v_y$ are provided as a finite time data set
and are interpolated using bicubic interpolation in space and 3rd order Lagrange polynomials in time.
This method
has been reported to be good enough for integrating trajectories in \cite{cf}.
We will focus our analysis in the time interval $[0, 900]$ in the area marked by a rectangle in Fig. \ref{fig:sfqgn}
for which \cite{nlpg} reports the computations of
several DHTs. In \cite{nlpg}  distinguished trajectories are computed by means of an iterative algorithm
which is initialized  on a hyperbolic instantaneous stagnation point (ISP). In particular two
 paths of such  ISPs are chosen in the Northern gyre in the time intervals $[0, 339]$ and $[446, 880]$.
From each  of these paths, a DHT is computed which is
in the same geographical area although its coordinates are determined for a  different time range.
In Fig. \ref{fig:dhtqgnab} we show the $x$ and $y$ evolutions for these trajectories. These coordinates have been computed
with a different algorithm  to that proposed in \cite{nlpg}. Instead each corresponds to a trajectory
which is in the intersection of  a piece of a stable manifold and a piece of an unstable manifold which are
evolved in  backwards and forwards time respectively.
In this procedure, in order to avoid the numerous intersections
between stable and unstable manifolds, which make difficult the tracking of the trajectory which is distinguished,
 manifolds are trimmed at each time step following the ideas  in \cite{physrep} where
a method is described to compute a piece of single
branch of  the stable or of the unstable manifold. This method  takes advantage of the fact that a
DHT must be in the intersection of both manifolds at all times, as it is a trajectory, however
does not improve the method
explained in  \cite{nlpg}  in the sense that it does not
allow either to extend the computation of  the DHT beyond the time interval in which the ISP exists.
Many questions have been raised for these trajectories as has been discussed in \cite{physrep,nlpg}.
For instance as
they have been computed only in finite time intervals on which the ISP exists, one can ask how to pursue
its computation beyond that interval.
Another open issue in \cite{nlpg} concerns deciding if  the two DHT  in Fig. \ref{fig:dhtqgnab}
computed at different times
are part of the same  trajectory.
In \cite{physrep}, the question is raised of whether it can happen that a DHT ceases to be distinguished or hyperbolic. In this section we apply our algorithm to compute limit coordinates
and verify whether trajectories in Fig. \ref{fig:dhtqgnab} are distinguished or not following our definition
\ref{def:dtdef}. Also we will describe how this definition helps address the questions raised in  \cite{physrep,nlpg}.
We have applied our algorithm to compute limit coordinates in the domain in which the DHT shown
in Fig. \ref{fig:dhtqgnab}a) exists. In particular we have applied it with the input:
$\mathcal{D}=(55, 75)\times (1325, 1375)$ km$^2$,  $t_0=120$, $t_N=300$, $\Delta t=5$ days, $\tau_0=2$ days, $\Delta \tau=5$ days, and
$\delta=10^{-3}$ km. The time step of the Runge-Kutta method is $0.1$ days.
Fig. \ref{fig:dhti}a) indicates with a solid line
the projection onto the $x-y$ plane of the trajectory depicted in figure \ref{fig:dhtqgnab}a)
in the interval $(120, 300)$, and with circles the path of limit coordinates. Fig. \ref{fig:dhti}b)
shows the evolution of the distances between these trajectories. This confirms that the trajectory displayed in Fig.
\ref{fig:dhtqgnab}a) is also distinguished in the sense of definition \ref{def:dtdef}
in the time interval $[120, 330]$ with accuracy $\epsilon=8 \cdot 10^{-1}$ km.
Thus in this time interval, limit coordinates give
a method for computing DT  different from those proposed in \cite{nlpg,kayo}. Circles in Fig. \ref{fig:dhtiellhyp}
show the location  versus time
of the $x$ limit coordinates computed with our algorithm. The solid line represents
a trajectory obtained after
integrating with  a 5th order Runge-Kutta method  forwards and backwards in time the initial condition of the circle
at day 285. The dashed line represents the same, but with the initial condition slightly perturbed.
It is evident that in both cases the trajectories are aligned with the path of limit coordinates.
The distinguished trajectory is highly hyperbolic backwards in time as in that
direction  a small perturbation
 amplifies greatly, while it does not do so forwards in time, suggesting that  it has a non-hyperbolic type of stability
in that direction (see comments to Figs. \ref{fig:duffingtray} and \ref{fig:duffingetray}).

Beyond day 300 it is possible to continue the path of limit coordinates. Figure
\ref{fig:dhtifin1} shows a diagram at day 330;  there  is shown the convergence of the
$x$ component of the minimum of $M$ versus $\tau$. This type of convergent diagram is not found
in this neighbourhood for day 337. On the other hand,  although it is possible to continue the path of limit coordinates
beyond day 300, Fig. \ref{fig:dhtifin2} proves that this path is not a trajectory. There can be seen
the existence of different trajectories crossing the path, confirming that it is not a trajectory
as otherwise it would violate the uniqueness of the solution. Therefore, following our construction
it is possible to say that beyond day 300 the trajectory is no longer distinguished.

Fig. \ref{fig:dhtd} confirms that the trajectory
in Fig. \ref{fig:dhtqgnab}b) is also distinguished in the sense of definition \ref{def:dtdef}
in the time interval $(470, 860)$ with accuracy $\epsilon=3$ km.
In particular to compute the path in Fig. \ref{fig:dhtd} we have applied the  algorithm of section
3 with the input:
$\mathcal{D}=(50, 65)\times (1255, 1270)$ km$^2$,  $t_0=470$, $t_N=860$, $\Delta t=5$ days, $\tau_0=2$ days, $\Delta \tau=7$ days, and
$\delta=10^{-3}$ km. The Runge-Kutta time step is $0.1$ days. In the time interval from day 600 to day 650
some of the input parameters were
 modified as follows: $\mathcal{D}=(73.5, 75.5)\times (1384, 1392)$ km$^2$, $\tau_0=40$ days, and $\Delta t=1$ day.
This  was due to the presence of nearby elliptic type minima in the function $M$, that made
difficult  tracking  the path of the limit coordinates with the previous input.

Finally, we discuss the existence of non-hyperbolic distinguished trajectories in this data set.
The presence of this type of trajectories has not been addressed before, and  we do not have any benchmark solution.
We have looked for this type of trajectory
in areas of the flow where Eulerian eddies seemed to persist for long times. Figure \ref{fig:Me}
represents the function $M$ at day 370  for
$\tau=150$ and $\tau=250$. In these figures there can be seen the structure of an eddie
 at the centre even for rather long $\tau$-values, however    figure \ref{fig:xytau370e} does not confirm the convergence
of the minimum of $M$ towards a constant value. On the other hand, the slow convergence
in diagram \ref{fig:duffingexy} towards the non-hyperbolic trajectory, already suggested that long time intervals
were required for that purpose, and those intervals might be difficult to find in realistic flows such
like the one analysed here, in which one is provided just with a finite time datat set.

\section{Conclusions}

In this paper we have proposed a new definition of distinguished trajectory that
attempts to extend the concept of fixed point and periodic orbit to aperiodic dynamical systems.
The concept of fixed point is trivially contained in the definition. Regarding other
 especially useful trajectories  in dynamical systems, for instance periodic orbits, we have not proved
that they fall within
the definition in a general way, but we have  numerically verified it for  selected 2D and 3D  examples.
The definition can be  implemented numerically
and  the procedure consists in determining a path of {\it limit coordinates}.
We have analysed exact examples for the Duffing equation with known distinguished trajectories,
both periodic an aperiodic, and we have found that the path of
  limit coordinates coincides, to within  numerical accuracy, with  the distinguished trajectories
 and therefore those trajectories are identified also as {\it distinguished } in the framework
of our definition. Our definition is novel with respect to previous works dealing with distinguished trajectories,
because it is applicable to non-hyperbolic trajectories. In particular we have studied
a periodic orbit of the Duffing equation with non-hyperbolic  stability and is also recognised as
{\it distinguished} by our definition. In this case the function $M$ from which
the {\it limit coordinates} are computed seems to be a suggestive tool for characterising Lagrangian eddies.
We have tested our definition in the context of realistic aperiodic flows where
distinguished hyperbolic trajectories had been found \cite{physrep,nlpg}. Again we have identified these trajectories
by paths of limit coordinates in certain time intervals. Beyond these time intervals the trajectories
are no longer {\it distinguished} according to our definition. Thus in the context of the
definitions provided in this paper, the property of a trajectory of being {\it distinguished} may be lost in time.
Also we have found evidence  that the hyperbolicity of
these trajectories is not  constant in time. These two statements
provide answers to the open questions mentioned in the text that have been addressed in \cite{physrep,nlpg}.

\section*{Acknowledgements}
The authors have been supported by CSIC grants PI-200650I224 and OCEANTECH (No. PIF06-059), Consolider I-MATH (C3-0104)
and the Comunidad de Madrid project SIMUMAT S-0505-ESP-0158.
 It is a pleasure to acknowledge many
useful conversations with Steve Wiggins, Des Small, Peter Haynes,  Emilio Hern\'andez-Garc\'{\i}a, Crist\'obal L\'opez, Antonio Turiel,
Emilio Garc\'{\i}a-Ladona, Michal Branicki, Wenbo Tang,  J.J. L. Vel\'azquez on numerous issues  related to this project.
We also thank very valuable suggestions from Daniel Fox, Andrew Thompson,
and Jodie Holdway.
The computational part of this work was done using the CESGA computers SVGD and FINIS TERRAE and using the
SIMUMAT-CSIC cluster ODISEA.


text:
\begin{document}

\maketitle
\begin{abstract}
Linear Feedback Shift Registers (LFRS) are tools commonly used in cryptography in many different context, for example as pseudo-random numbers generators. In this paper we characterize LFRS with certain symmetry properties. Related to this question we also classify polynomials $f$ of degree $n$ satisfying the property that if $\alpha$ is a root of $f$ then $f(\alpha^n)=0$. The classification heavily depends on the choice of the fields of coefficients of the polynomial; we consider the cases $K=\F_p$ and $K=\Q$.
\end{abstract}

\section{Introduction}

The motivation of this paper comes from an exercise in a written exam of Cryptography about the bit stream $(s_j)$ ($j=0,1, \cdots$) generated by a (Fibonacci $n$-bit) Linear Feedback Shift Register (LFSR) see e.g. \cite[page 374]{Schneier15}:

\begin{figure}[h!]
\begin{center}
\includegraphics[height=3cm,width=10cm]{LFSRm.pdf}
\end{center}
\caption{a Fibonacci LFSR}
\end{figure}

A linear feedback shift register (LFSR) is a shift register whose input bit is a linear function of its previous state.
The bit stream
\[
\cdots s_n \, s_{n-1} \, \cdots \, s_2 \,  s_1 \,  s_0
\]
is generated recursively as in the Fibonacci sequence; namely, given the initialization state of the register $s_{n-1}, \cdots, s_0 \in \F_2$, for every $j \geq n$
 the bit $s_j$ is computed recursively as
\[
s_j = s_{j-1} p_{n-1} + s_{j-2} p_{n-2} + \cdots + s_{j-n} p_0  \, \, (\mathrm{mod} \, \, 2).
\]
The polynomial $\chi(x) = x^n + p_{n-1}x^{n-1} + \cdots + p_1 x + p_0$ characterizes the LFSR.

Because the operation of the register is deterministic, the stream of values produced by the register is completely determined by its current (or previous) state. Likewise, because the register has a finite number of possible states, it must eventually enter a repeating cycle. However, a LFSR with a well-chosen feedback function can produce a sequence of bits which has a very long cycle and so appears random. Applications of LFSRs include generating pseudo-random numbers, pseudo-noise sequences, fast digital counters, and whitening sequences.

\medskip

By using a row $\mathrm{w}=[s_{n-1} s_{n-2} \cdots s_1 s_0 ]$ to describe the state of the register, we have that
it changes according to a right multiplication $ \mathrm{w} \cdot \mathrm{L}$, where $\mathrm{L}$ is the $n \times n$ matrix
\begin{equation}\label{LFSRm}
\mathrm{L} =
\begin{bmatrix}
{ p_{n-1}}  & {1} & 0 & 0 & \cdots & 0 \\
{ p_{n-2}}  & 0 & {1} & 0 & \cdots & 0 \\
{ p_{n-3}}  & 0 & 0 & {1} & \cdots & 0 \\
\vdots  &\vdots & \vdots & \vdots & \ddots & \vdots\\
{ p_1 } & 0 & 0 & 0 & \cdots & {1} \\
{ p_0 } & 0 & 0 & 0 & \cdots & 0 \\
\end{bmatrix}
\end{equation}

Notice that $\chi(x)$ is exactly the characteristic polynomial of $\mathrm{L}$.

\medskip

The above mentioned exercise asked to run the $2$-bit LFSR with polynomial $\chi(x)=x^2 + x + 1$ and to compute the first $6$ bits
$s_5,s_4,s_3,s_2,s_1,s_0$ given $s_1 = 0 , s_0 = 1$.
What grabbed our attention and motivated this paper is that a student constructed the correct bit stream but in a wrong way. Namely, he constructed a bit stream $(r_j)$ by computing \[ [r_j \, r_{j+1} ] = [r_{j-2} \,  r_{j-1}] \cdot  \mathrm{L} \, .\]
It turns out that the bit streams $(r_j)$ and $(s_j)$ are the same because, in this specific case, the matrix associated to the LFSR satisfies
\begin{equation}{\label{2-bit}}
 \mathrm{L} = \tau \cdot \mathrm{L}^2 \cdot \tau,
\end{equation}
where $\tau$ is the $2 \times 2$ reflection matrix $\tau = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$.

 It is then natural to ask, more in general, about the classification of $n$-bit LFSR such that
 \begin{equation}{\label{n-bit}} \mathrm{L} = \tau \cdot \mathrm{L}^n \cdot \tau, \end{equation}
where $\tau = \begin{bmatrix} 0 & \cdots & 1 \\ 0 & \iddots & 0 \\ 1 & \cdots & 0 \end{bmatrix}$ is the $n \times n$ reflection matrix. We will prove the following result:

\begin{theorem}\label{thm:LFSRthm} Let $\mathrm{L}$ be a $n \times n$ matrix with entries in $\F_2$ of the form (\ref{LFSRm}) and let $\chi(x)$ be its characteristic polynomial. Then, $\mathrm{L}$ satisfies equation (\ref{n-bit}) if and only if $\chi(x)= \sum_{j=0}^{n} x^j$.
\end{theorem}

More in general, \eqref{n-bit} implies that the matrices $\mathrm{L}$ and $\mathrm{L}^n$ are similar, and so they have the same eigenvalues. In particular, this implies that for every root $\alpha \in \overline{\F}_2$ of $\chi(x)$, then $\chi(\alpha^n)=0$. This led to the following natural question: given a field $K$, classify all the polynomials $f(x)\in K[x]$ of degree $n$ satisfying the property
\begin{equation} \label{main_property}
\mbox{for every } \alpha \in \overline{K} \mbox{ such that } f(\alpha)=0, \mbox{ then } f(\alpha^n)=0.
\end{equation}
Notice that the polynomial $\chi(x)=\sum_{j=0}^n x^j$, arising from Theorem \ref{thm:LFSRthm}, satisfies this property.

Trivially, every linear polynomial satisfies \eqref{main_property}, so we will always assume without loss of generality that $\deg f \ge 2$.
If $f$ satisfies the property, then for every $i\ge 0$, $\alpha^{n^i}$ is a root of $f$. As $f$ has at most $n$ distinct roots, this implies that there exist $0\le h < k \le n$ such that $\alpha^{n^k}= \alpha^{n^h}$. If we rewrite this equality as $\alpha^{n^h}(\alpha^{{n^h}(n^{k-h}-1)}-1)=0$, then we have that either $\alpha=0$ or $\alpha$ is a root of unity. This property is independent of the nature of the field of coefficients $K$, while the characterization heavily depends on it.

In this paper we will give the full classification in the case  $K=\Q$ (Section \ref{sec:rational}) and for irreducible polynomials in the case $K=\F_p$ (Section \ref{sec:finite_fields}). In both cases, these polynomials are strictly connected to the cyclotomic polynomials, and in particular it turns out their degrees are connected with other classical problems in number theory, e.g. the classification of the numbers which are coprime with their Euler totient function, which arise in several different contexts, as explained in Section \ref{sec:Michon_numbers}.

\subsection*{Acknowledgements}
Both the authors are members of the INdAM group GNSAGA, of CrypTO (the group of Cryptography and Number Theory of Politecnico di Torino), and of DISMA, Dipartimento di Eccellenza MIUR 2018-2022.

\section{Proof of Theorem \ref{thm:LFSRthm}}

\noindent In this section we prove Theorem \ref{thm:LFSRthm}.

\medskip

\noindent First, assume that $\mathrm{L} = \tau \cdot \mathrm{L}^n \cdot \tau$; then  $\tau \cdot \mathrm{L} = \mathrm{L}^n \cdot \tau$.
A direct computation shows that:
\[
\tau \cdot \mathrm{L} =
\begin{bmatrix}
{ p_{0}}     & 0 & * & * & \cdots & * \\
{ p_{1}}    & 0 & * & * & \cdots & * \\
{ p_{2}}    & 0 & * & * & \cdots & * \\
\vdots      &\vdots & \vdots & \vdots & \ddots & \vdots\\
{ p_{n-2} } & 0 & * & * & \cdots & * \\
{ p_{n-1} } & 1 & * & * & \cdots & * \\
\end{bmatrix},
\]
and
\[
\mathrm{L}^n \cdot \tau =
\begin{bmatrix}
{ p_{n-1}} & p_{n-2} + p_{n-1}^2 & * & * & \cdots & * \\
{ p_{n-2}} & p_{n-3} + p_{n-2} p_{n-1} & * & * & \cdots & * \\
{ p_{n-3}} & p_{n-4} + p_{n-3} p_{n-1}  & * & * & \cdots & * \\
\vdots     &\vdots & \vdots & \vdots & \ddots & \vdots\\
{ p_{1} }  & p_0 + p_1 p_{n-1} & * & * & \cdots & * \\
{ p_{0} }  & p_0 p_{n-1} & * & * & \cdots & * \\
\end{bmatrix},
\]
which holds in any characteristic.
From these formulas, as the coefficients of the characteristic polynomial lie in $\F_2$, we directly get that $p_0=p_1=\cdots=p_{n-1}=1$, i.e. $\chi(x)=\sum_{j=0}^n x^j$.
\medskip

\noindent Conversely, assume that $\chi(x)=\sum_{j=0}^n x^j$ i.e. \[ \mathrm{L} =
\begin{bmatrix}
{ 1}  & {1} & 0 & 0 & \cdots & 0 \\
{ 1}  & 0 & {1} & 0 & \cdots & 0 \\
{ 1}  & 0 & 0 & {1} & \cdots & 0 \\
\vdots  &\vdots & \vdots & \vdots & \ddots & \vdots\\
{ 1 } & 0 & 0 & 0 & \cdots & {1} \\
{ 1} & 0 & 0 & 0 & \cdots & 0 \\
\end{bmatrix} \, .\]
Let us write $\mathrm{L} = [E_0 E_1 E_2 \cdots E_{n-1}]$, where $E_j$ denote the $(j+1)^{\rm{th}}$ column of $\mathrm{L}$. Notice that, for every matrix $\mathrm{M} = [M_1 M_2 \cdots M_n]$, then the following holds:
\[ \mathrm{M} \cdot \mathrm{L} = \left [ \left ( \sum_j M_j \right ) \, M_1 \, M_2 \, \cdots \, M_{n-1} \right ] \, . \]
Using this property with $M=L$, it follows that
\[\mathrm{L}^2 = [E_n \, E_0 \, E_1 \, \cdots \, E_{n-2}] \]
where $E_n = (\sum_{j=0}^{n-1} E_j)$.
Since we are working in $\F_2$, a simple inductive argument implies that, for every $3 \le k \le n$,
\[
\mathrm{L}^k = [E_{n-k+2} \, E_{n-k+3} \, \cdots \, E_n \, E_0 \, \cdots \, E_{n-k-1}];
\]
in particular
\[
 \mathrm{L}^n = [E_2 E_3 \cdots E_{n-1} E_n E_0] \, .
\]
Now it is straightforward to check that $ \tau \cdot \mathrm{L} = \mathrm{L}^n \cdot \tau$, as wanted. \qed

\begin{remark}
We point out that the statement of Theorem \ref{thm:LFSRthm} heavily depends on the fact that the field of the coefficients has characteristic $2$. More in general, it can be proved in the same way that, if char$(K)\neq 2$, a matrix $L$ of the form \eqref{LFSRm} satisfies $\mathrm{L}=\tau \cdot \mathrm{L}^n \cdot \tau$ if and only if
\begin{equation*}
\mathrm{L} =
\begin{bmatrix}
{-1}  & {1} & 0 & 0 & \cdots & 0 \\
{-1}  & 0 & {1} & 0 & \cdots & 0 \\
{-1}  & 0 & 0 & {1} & \cdots & 0 \\
\vdots  &\vdots & \vdots & \vdots & \ddots & \vdots\\
{-1} & 0 & 0 & 0 & \cdots & {1} \\
{-1} & 0 & 0 & 0 & \cdots & 0 \\
\end{bmatrix}
\end{equation*}
Notice that, if  char$(K)\neq 2$, then a matrix $\mathrm{L}$ of the form \eqref{LFSRm} has characteristic polynomial equal to $\chi(x)=x^n-p_{n-1}x^{n-1}-\cdots -p_0$. This implies that the previous matrix has characteristic polynomial $\chi(x)=\sum_{i=0}^n x^i$, which is a product of cyclotomic polynomials.

\end{remark}

\section{Cyclotomic polynomials} \label{sec:cyclotomic}

In this section we recall the definition and some basic properties of cyclotomic polynomials which will be used later. For more references, see \cite{Lang} or \cite{finite_fields}.

Given an integer $n\ge 0$, we let $\Phi_n(x)$ denote the $n^{th}$-cyclotomic polynomial, defined by
\[
\Phi_n(x):=\prod_{\substack{ 1\le k \le n \\ (k,n)=1} }\left (x- e^{\frac{2\pi i k}{n}} \right ).
\]
Clearly, we have that $x^n -1 = \prod_{d \mid n} \Phi_d (x)$ for every $n \ge 1$, and the M\"obius inversion formula gives
\[
\Phi_n(x)=\prod_{d \mid n} \left ( x^d-1 \right )^{\mu(n/d)}= \prod_{d \mid n} \left ( x^{n/d}-1 \right )^{\mu(d)},
\]
where $\mu(d)$ denotes the M\"obius function, i.e.
\[
\mu(d):= \begin{cases} 1 \ \ \quad \mbox{if } d \mbox{ is the product of a even number of distinct prime factors;}		\\
0 \ \ \quad \mbox{if } d \mbox{ is not squarefree;}\\
-1 \quad \mbox{if } d \mbox{ is the product of a odd number of distinct prime factors.}
\end{cases}
\]
It can be proved that, for every $n \ge 0$, $\Phi_n$ is a monic polynomial with integer coefficients and that $(\Phi_m(x), \Phi_n(x))=1$ for every $m<n$. The degree of $\Phi_n$ is clearly equal to $\varphi(n)$, which denotes the Euler totient function, i.e.
\[
\varphi(n)= n \prod_{p \mid n} \left ( 1- \frac{1}{p} \right ).
\]

\begin{example}
If $n=p$ is a prime, then $\Phi_p(x)=x^{p-1}+ x^{p-2} + \cdots + 1$.
\end{example}

Let $p$ be a prime and let $\F_q$ denote the finite field with $q=p^f$ elements for some $f\ge 1$. The cyclotomic polynomials are irreducible over the field of rational numbers, but this is not the case over finite fields; indeed, the polynomial $\Phi_n$ is irreducible over $\F_q$ if and only if $n=2, 4, r^k$ or $2r^k$, where $r$ is an odd prime, $k$ is a positive integer and $q$ is a generator of $(\Z/n \Z)^{\times}$.\footnote{In particular this implies that if $(\Z/n\Z)^{\times}$ is not cyclic (i.e. unless $n$ is an odd prime power, twice an odd prime power, or $n=2$ or $4$), then $\Phi_n$ is a polynomial which is reducible modulo every prime $p$ but is irreducible over $\Q$.}

If $(n,q)=1$, then $\Phi_n$ can be factorized into $\varphi(n)/m$ distinct irreducible polynomials of the same degree $m$ over $\F_q$, where $m$ is the multiplicative order of $q$ modulo $n$, i.e. the least positive integer such that $q^m \equiv 1 \mod n$. For the rest of the paper we will denote this order by $\ord_n(q)$.
 If $(n,q)\neq 1$ we can write $n=p^a n'$ with $(p,n')=1$; then, we have that $\Phi_{n}(x)=(\Phi_{n'}(x))^{p^a}$, so it is enough to study the factorization for cyclotomic polynomials of roots of unity of order coprime with the characteristic of the field.

\begin{example} \label{ese:1}
Let us consider the cyclotomic polynomial of order $15$ over $\F_2$; we have that $\varphi(15)=8$ and $\ord_{15}(2)=4$, hence $\Phi_{15}$ factorises into $2$ factors of degree $4$, i.e.
\[
\Phi_{15}(x)=(x^4+x^3+1)(x^4+x+1).
\]
If we take the cyclotomic polynomial of order $8$ over $\F_2$, then we have that
\[
\Phi_8(x)=(x+1)^4.
\]

\end{example}

\section{Classifying the polynomials with rational coefficients satisfying property \eqref{main_property}}

Let $f(x) \in K[x]$ be a nonzero polynomial of degree $n$ with coefficients in a field $K$; we are interested in characterizing the polynomials with the following property:
\begin{equation*}
\mbox{for every } \alpha \in \overline K \mbox{ such that } f(\alpha)=0, \mbox{ then } f(\alpha^n)=0.
\end{equation*}

We will be interested in two main cases, namely the case $K=\F_p$ and the case $K=\Q$. For every choice of the field $K$, we will first look at the irreducible case and then at the general case. Notice that, eventually dividing by the leading coefficient of $f$, it is enough to consider $f$ monic.

\subsection{The case of irreducible polynomials over $\Q$}

We saw at the end of the introduction that trivially all the polynomials of degree $1$ satisfy the condition, so without loss of generality we can assume $\deg f \ge 2$. In this case, if $f$ is a polynomial of degree $\ge 2$ having this property, then its roots are either $0$ or roots of unity. If $f$ is irreducible and monic, this means that $f(x)$ is the minimal polynomial of a (primitive) root of unity, i.e. a cyclotomic polynomial.

We can prove the following proposition, which gives the characterization of the irreducible polynomials which satisfy \eqref{main_property}.

\begin{proposition} \label{prop:irreducible}
Let $f$ be a monic, irreducible polynomial of degree $n\ge 2$; then, $f$ satisfies \eqref{main_property} if and only if $f(x)=\Phi_d(x)$ for some $d\ge 3$ with $(n,d)=1$.
\end{proposition}

\begin{proof}
If $f$ satisfies the property \eqref{main_property}, then for every root $\alpha$ of $f$, either $\alpha=0$ or $\alpha$ is a root of unity.

As we are assuming $f$ irreducible of degree $\ge 2$, then $f(x)$ is the minimal polynomial of a root of unity, i.e. a cyclotomic polynomial. Moreover, if the order of $\alpha$ is $d$, then $f(x)=\Phi_d(x)$ and $\deg f=\varphi(d)$, so we have to take $d\ge 3$ because we are assuming that $d\ge 2$. As we know that all the roots of $\Phi_d$ are primitive $d^{th}$-roots of unity, then $\alpha^n$ has to be a primitive $d^{th}$ root of unity, which implies that $(n,d)=1$, proving the first implication.
\medskip

We want now to prove the converse. Let us assume then that $f(x)=\Phi_d(x)$ is the $d^{th}$-cyclotomic polynomial with $d\ge 3$ and that the degree of $f$ is coprime with $d$. First, as $d \ge 3$, then $\deg f=\phi(d)$ as wanted. Moreover, if $\alpha$ is a root of $f$, then it is a primitive $d^{th}$-root of unity, and as $(n,d)=1$, we have that $\alpha^n$ is also a primitive $d^{th}$-root of unity, and so $f(\alpha^n)=0$, proving that $f$ satisfies the condition \eqref{main_property}, as wanted.
\end{proof}

\begin{example}\label{ese_sec4} \
\begin{itemize}
  \item If $n=p-1$ where $p$ is a prime, then $\Phi_p(x)=x^{p-1}+ \cdots + 1$ satisfies the condition \eqref{main_property};
  \item Let $\deg f=60$; there are exactly two polynomials of degree $60$ which are cyclotomic polynomials and satisfy the condition \eqref{main_property}, i.e. $\Phi_{61}(x)$ and $\Phi_{77}(x)$ (and $60$ is the smallest degree with this property).
 \end{itemize}
\end{example}

We just proved that the irreducible polynomials which satify the condition \eqref{main_property} are exactly the polynomials of degree $1$ and the cyclotomic polynomials $\Phi_d(x)$ with $(\deg f, d)=1$. In particular, since $\deg \Phi_d(x)=\varphi(d)$, where $\varphi$ is the Euler totient function, we are asking that $(d, \varphi(d))=1$. This also implies that if $\deg f >1$, then it has to be odd and squarefree (see Proposition \ref{prop:cyclic_char}).

The integers $d$ such that $(d, \varphi(d))=1$ have been deeply studied in the literature and appear in many different contexts. We are going to describe some properties of these numbers Section \ref{sec:Michon_numbers}.

\subsection{The general case over $\Q$} \label{sec:rational}

We want now to analyse the general case and characterise the polynomials $f$ (not necessarily irreducible) such that $f$ satisfies the property \eqref{main_property}.

As seen before, if $\alpha$ is a root of $f$ then either $\alpha$ is zero or $\alpha$ is a root of unity, hence either $x \mid f$ or, if $\alpha$ is a primitive $k^{th}$ root of unity, $\Phi_k \mid f$.

We give the following general characterization for the polynomials which satisfy property \eqref{main_property}:

\begin{theorem} \label{thm:general}
Let $f(x)\in \Q[x]$ be a monic polynomial of degree $n\ge 2$. Then, $f$ satisfies the property \eqref{main_property} if and only if
\begin{equation} \label{eq:shape}
f(x)=x^a \prod_{\substack{ 1\le h_1 < \cdots < h_r \le n \\ (n,h_i) = 1}} \Phi_{h_i}(x)^{b_i} \prod_{\substack{ 1\le k_1 < \cdots < k_s \le n \\ (n,k_j) \neq 1}}
 \Phi_{k_j}(x)^{c_j} \prod_{t=1}^{m_j} \Phi_{\frac{k_j}{(n^t, k_j)}}(x)^{d_{t,j}},
\end{equation}
where $h_i, k_j$ are positive integers, $a,b_i, c_j, d_j $ are non-negative integers, $m_j$ is the biggest integer $\le n$ such that $\frac{k_j}{(n^t, k_j)}$ is not coprime with $n$ and, if $c_j\neq 0$, then also $d_{t,j} \neq 0$ for every $t=1, \ldots, m_j$. Moreover
\begin{equation} \label{eq:degree}
 n= a+\sum_{i=1}^r b_i \varphi(h_i) + \sum_{j=1}^s \left ( c_j\varphi(k_j) + \sum_{t=1}^{m_j} d_{t,j} \varphi \left (\frac{k_j}{(n^{t},k_j)}\right ) \right ).
 \end{equation}
\end{theorem}

\begin{remark}
We notice that if $f$ is irreducible, $f$ has only one irreducible factors, so it is equal to some $\Phi_d(x)$ with $d\ge 3$, $n=\varphi(d)$ and $(n,d)=1$, as proved in Proposition \ref{prop:irreducible}.
\end{remark}

\begin{proof}
Assume first that $f$ is a monic polynomial of degree $n\ge 2$ satisfying the property \eqref{main_property}; as its roots are either $\alpha=0$ or $\alpha$ a root of unity then, then the irreducible factors of $f$ are either $x$ or cyclotomic polynomials. Assume that $\alpha \in \Qbar$ is a root of $f$ which is a primitive root of unity and denote by $k$ its order. As $\Phi_k(x)$ is the minimal polynomial of $\alpha$, then $\Phi_k \mid f$. If $(n,k)=1$, then $\alpha^n$ is again a primitive $k$-th root of unity, hence its minimal polynomial is still $\Phi_k$. Assume that $(n,k)\neq 1$ and let $m$ be the biggest integer $\le n$ such that $k/(n^m,k)$ is not coprime with $k$. Then, for every $t=1, \ldots, m$, we have that $\alpha^{n^t}$ is a primitive $k/(n^t,k)$-th root of unity, hence $\Phi_{k/(n^t,k)}$ must divide $f$. On the other hand, if $k/(n^t,k)$ is coprime with $n$, then $(n^t,k)=(n^m,k)$ and so $\alpha^{n^t}$ is a $k/(n^{m},k)$-th primitive root of unity for every $u=m+1, \ldots, n$, hence its minimal polynomial is again $\Phi_{k/(n^{m},k)}$. This implies that $f$ has to be of the form \eqref{eq:shape}. Moreover, the relation \eqref{eq:degree} comes directly by computing the degree of the product, recalling that $\deg \Phi_k(x)=\varphi(k)$. \\

We want to prove the converse; assume that $f$ is a polynomial of the shape \eqref{eq:shape}; then, the degree $n$ of $f$ satisfies \eqref{eq:degree}. Consider now a root $\alpha$ of $f$; then, $\alpha$ is a root of one of the irreducible factors of $f$. From \eqref{eq:shape}, then either $\alpha$ is equal to $0$, or $\alpha$ is a root of a cyclotomic polynomial, i.e. it is a root of unity.

If $\alpha=0$, then $\alpha^n=0$ and so $f(\alpha^n)=0$ as wanted. Assume now that $\alpha$ is a root of unity and denote by $k$ its order. As remarked before, if $\alpha$ has order $k$ then $\alpha^{n^t}$ has order $k/(n^t,k)$ for every $t=1, \ldots, n$. This implies that, if $(n,k)=1$ then $\alpha^{n^t}$ is a primitive $k$-root of unity, so $\Phi_k(\alpha^{n^t})=0=f(\alpha^{n^t})$. Assume now that $(n,k)\neq 1$; then, for every $t$ also $\Phi_{\frac{k}{(n^t,k)}} \mid f$ and so $f(\alpha^{n^t})=\Phi_{\frac{k}{(n^t,k)}}(\alpha^{n^t})=0$, as wanted.
\end{proof}

\begin{example}
For every $n \ge 1$, the polynomial $x^n-1$ satisfies the property \eqref{main_property}. Notice that
$$ x^n-1= \prod_{d \mid n} \Phi_d(x), $$
so $x^n-1$ has exactly the shape \eqref{eq:shape}.
\end{example}

\begin{example}
Let us list all the polynomials of degree $6$ which satisfy the property \eqref{main_property}.
By Theorem \ref{thm:general}, the possible factors of $f$ are either $x$ or cyclotomic polynomials of degree $\le 6$. We use a result of Gupta \cite{Gupta}, which asserts that if $n \in \varphi^{-1}(m)$, then $n < m < A(m)$ with
$$ A(m)= m \prod_{p-1 \mid m} \frac{p}{p-1}. $$
Using this, we have that $A(6)=21$, so we have to check the cyclotomic factors up to $\Phi_{21}$. An easy calculation gives that the cyclotomic polynomials with degrees equal to $1$ are $\Phi_1$ and $\Phi_2$, the ones of degree $2$ are $\Phi_i$ with $i=2,3,6$, the ones of degree $4$ are $\Phi_i$ with $i=5,8,9,12$ and the ones of degree $6$ are $\Phi_i$ with $i=7,14,18$. In order to classify the possible polynomials satisfying \eqref{main_property}, we have to take into account the fact that, if $\Phi_k \mid f$ for some $k$, then also $\Phi_{\frac{k}{(6,k)}} \mid f$.
The polynomials $f$ of degree $6$ that satisfy the property \eqref{main_property} have one of the following shape:
\begin{itemize}
   \item $f(x)=x^a\Phi_1(x)^b$ with $a+b=6$;
   \item $f(x)=x^a\Phi_1(x)^b\Phi_2(x)^c$ with $a+b+c=6$ and $b,c \neq 0$;
   \item $f(x)=x^a\Phi_1(x)^b\Phi_2(x)^c\Phi_3(x)^d$ with $a+b+c+2d=6$ and $b,d \neq 0$;
   \item $f(x)=x^a\Phi_1(x)^b\Phi_2(x)^c\Phi_3(x)^d\Phi_4(x)^{e}$ with $a+b+c+2d+2e=6$ and $b,d,e \neq 0$;
   \item $f(x)=x^a\Phi_1(x)^b\Phi_2(x)^c \Phi_5(x)$ with $a+b+c=2$;
   \item $f(x)=x^a \Phi_1(x)^b\Phi_2(x)^c\Phi_3(x)^d\Phi_4(x)^{e} \Phi_6(x)^f$ with $a+b+c+2(d+e+f)=6$ and $b,f \neq 0$;
   \item $f(x)=\Phi_7(x)$ (which is the only irreducible $f$);
   \item $f(x)=\Phi_1(x) \Phi_2(x) \Phi_{12}(x)$.
\end{itemize}
\end{example}

\section{Numbers coprime with their Euler totient function} \label{sec:Michon_numbers}

The integers $d$ such that $(d, \varphi(d))=1$ have been studied in number theory and appear in many different contexts; for example, these are the numbers such that there is only one group of order $d$ (i.e. the cyclic one). For this reason the numbers which satisfy this property are usually called \textit{cyclic}.

If $n$ is a prime, then $\varphi(n)=n-1$ so $n$ is a cyclic number; this shows that cyclic numbers are infinite. In \cite{Erdos}, Erd\"os gave an asymptotic formula for the number of cyclic numbers.

We prove the following easy necessary condition:
\begin{proposition} \label{prop:cyclic_char}
If $d$ is a cyclic number then either $d=2$ or $d$ is odd and squarefree.
\end{proposition}

\begin{proof}
We first notice that, if $d> 2$, then $\varphi(d)$ is even, so the only even number $d$ which is cyclic is $d=2$. On the other hand, assume that $p^2 \mid d$; then $p \mid \varphi(d)$, and so $(d, \varphi(d))\neq 1$ which contradicts the property of being cyclic.
\end{proof}

Of course these conditions are not sufficient in general; in fact, if we take for example $d=21$, then $\varphi(21)=12$ which is not coprime with $21$.

We can however prove the following result:

\begin{proposition}
 If $d$ is an odd numbers which is the product of two consecutive prime numbers, then $d$ is cyclic.
\end{proposition}

\begin{proof}
Assume that $d=p_n p_{n+1}$, where $p_i$ denotes the $i^{th}$ prime number; then we can prove that $p_n \nmid (p_{n+1}-1)$. In fact, assume by contradiction that $p_n \mid (p_{n+1}-1)$; as $p_n \neq 2$, then $p_n \mid \frac{p_{n+1}-1}{2}$, but this is a contradiction since by Bertrand's postulate $p_n < p_{n+1} < 2p_n$.
\end{proof}

Cyclic numbers are also related to Carmichael numbers. We recall that Carmichael numbers \cite{Carmichael} are composite numbers $n$ which satisfies the modular arithmetic congruence relation:
\[
b^{n-1} \equiv 1 \pmod n
\]
for all integers $b$ which are relatively prime to $n$.  Carmichael numbers are also called Fermat pseudoprimes or absolute Fermat pseudoprimes. Indeed, Carmichael numbers pass a Fermat primality test with respect to every base b relatively prime to the number, even though it is not actually prime. This makes tests based on Fermat's Little Theorem less effective than strong probable prime tests such as the Baillie-PSW primality test and the Miller-Rabin primality test. Korselt \cite{Korselt} proved that a positive composite integer $n$ is a Carmichael number if and only if $n$ is square-free, and for all prime divisors $p$ of $n$ then $p-1 \mid n-1$. For Carmichael numbers the following proposition holds:

\begin{proposition}
Every divisor of a Carmichael number is odd and cyclic.
\end{proposition}

In the '80, Michon conjectured that the converse is also true, i.e. that every odd cyclic number has at least one Carmichael multiple. The conjecture has been verified by Crump and Michon  for all the numbers $< 10000$, but remains still open.

\begin{example}
Here we list the cyclic numbers $1< d < 100$ which are not prime and the relative Euler totient functions:
\begin{tabbing}
marrone \= bla \= bla \= bla \= bla \= bla \= bla \= bla \= bla \= bla \= bla \= bla \= bla \kill
$d$ \>  15   \> 33   \> 35  \> 51   \> 65   \> 69  \> 77  \> 85 \> 87   \> 91  \> 95 \\
$\varphi(d)$ \> 8 \> 20 \> 24 \> 32 \> 48 \> 44 \> 60 \> 64 \> 56  \> 72 \> 72 \\
\end{tabbing}
This shows for example that for the numbers $\le 100$ the only cyclic numbers which have the same Euler totient function are $61$ and $77$ (for which $\varphi(61)= \varphi(77)=60$, see Example \ref{ese_sec4}) and $73, 91$ and $95$ (for which $\varphi(73)= \varphi(91)= \varphi(95)=72$).
\end{example}

\section{The irreducible case over $\F_p$} \label{sec:finite_fields}

The case of finite fields is very different from the previous one. Also in this case, all linear polynomials satisfy property \ref{main_property}, so without loss of generality we will always assume that $\deg f \ge 2$. In this setting, it is always true that, if $\alpha\in \F_q$ for some $q=p^a$, then $\alpha$ has finite order. We will denote by $k$ its order, i.e. the minimal $k$ such that ${\alpha}^k =1$. Notice that $(k,p)=1$. If we consider the cyclotomic polynomial $\Phi_k$,  it is not true anymore that this is the minimal polynomial of $\alpha$ as cyclotomic polynomials are not always irreducible over $\F_p$. Indeed, $\Phi_k$ factorises into $\varphi(k)/\ord_k(p)$ irreducible polynomials of degree $n:=\ord_k(p)$, where $\ord_k(p)$ denotes the multiplicative order of $k$ modulo $p$. This means that, if we consider the extension $\F_p(\alpha)/\F_p$, then it has degree $n$. Moreover, as this is a finite field extension, it is Galois with cyclic Galois group of order $n$, and a generator of the Galois group is the Frobenius, i.e. the automorphism with sends $\alpha \mapsto \alpha^p$. This means that all the conjugates of $\alpha$ are exactly $\alpha, \alpha^p, \ldots \alpha^{p^{n-1}}$. We are ready to prove the following:

\begin{proposition}
Let $f$ be an irreducible polynomial of degree $n\ge 2$; then, $f$ satisfies property \eqref{main_property} if and only if $f$ is a factor of degree $n$ of a cyclotomic polynomial $\Phi_k$ with $n < \varphi(k)$ and $n$ is a power of $p$, or $f=\Phi_k$, ${\ord}_k(p)=\varphi(k)$ and $(k, \varphi(k))=1$.
\end{proposition}

\begin{proof}
First assume that $f$ is an irreducible polynomial satisfying property \eqref{main_property}; then, if $\alpha$ is a root of $f$, either $\alpha$ is zero or $\alpha$ is a root of unity. If $\alpha=0$, then $f(x)=x$, which we exclude as we are assuming $\deg f \ge 2$.  Hence $\alpha$ is a primitive $k$-th root of unity for some $k\in \N$. Notice that $(k,p)=1$ since we are in characteristic $p$ and so $a^p=a$ for every $a\in \overline{\F}_p$. As $f$ is irreducible, we have that $f$ divides the cyclotomic polynomial $\Phi_k$. Now, if $\Phi_k$ is irreducible over $\F_p$ (which, as seen in Section \ref{sec:cyclotomic}, happens if and only if $\ord_k(p)=\varphi(k)$), then $f=\Phi_k$. Moreover, by property \eqref{main_property} we have that $\alpha^{{\varphi(k)}^t}$ must be a primitive $k$th root of unity for every $t \ge 1$, which implies that $(k,\varphi(k))=1$ as wanted.

Assume now that $\ord_k(p)<\varphi(k)$; then, $\Phi_k$ factorises into $\varphi(k)/\ord_k(p)$ irreducible factors of degree $\ord_k(p)$ and $f$ will be equal to one of these factors. This implies that $[\F_p(\alpha):\F_p]=\ord_k(p)$ and, using the Frobenius, all the other roots of $f$ will be $\alpha^p, \ldots, \alpha^{p^{k-1}}$. Therefore, as by property \eqref{main_property} we have that $\alpha^{\ord_k(p)}$ is a root of $f$, this implies that $\ord_k(p)$ must be a power of $p$, as wanted.
\medskip

Let us prove the converse. Assume first that $f=\Phi_k$ with $\ord_k(p)=\varphi(k)$ and $(k, \varphi(k))=1$; then, it is easy to prove that $f$ satisfies property \eqref{main_property} since $\alpha^{\varphi(k)^t}$ is again a primitive $k$-th root of unity as $(k, \varphi(k))=1$. Let us finally consider the case in which $f$ is an irreducible factor of degree $n< \varphi(k)$ and $n$ is a power of $p$; as seen in Section \ref{sec:cyclotomic}, we have that $n=\ord_k(p)$. As by assumption $n$ is a power of $p$, this directly implies that $\alpha^{n^t}$ is a root of $f$ since and all the cojugates of $\alpha$ are exactly $\alpha^p, \ldots \alpha^{p^{n-1}}$, concluding the proof.
\end{proof}

The last proposition gives strong constraints on the type of polynomials which satisfies property \eqref{main_property}; in particular:
\begin{itemize}
\item if $\Phi_k$ is irreducible over $\F_p$, then $(k,p)=1$ and $(k,\varphi(k))=1$;
\item if $\Phi_k$ is not irreducible over $\F_p$ and $(k,p)=1$, then $\ord_k(p)$ must be a power of $p$.
\end{itemize}
In Section \ref{sec:cyclotomic} we recalled that $\Phi_k$ is irreducible over $\F_p$ if and only if $k=2,4, r^m$ or $2r^m$ with $r$ an odd prime, and the multiplicative order of $p$ modulo $k$ is maximal. By Proposition  \ref{prop:cyclic_char}, we have that if $(k,\varphi(k))=1$ then either $k=2$ or $k$ is odd and squarefree; combining these two conditions we have that $\Phi_k$ is irreducible over $\F_p$ with $(k,\varphi(k))=1$ if and only if either $k=2$ and $p \neq 2$ or $k$ is an odd prime different from $p$ and $\ord_k(p)=k-1$.\\

In the second case we have even a more restricted condition; indeed, $\ord_k(p)$ is a divisor of $\varphi(k)$. Let us write $k=\prod_{i=1}^m p_i^{a_i}$ where the $p_i$ are distinct primes and $a_i$ are positive integers; then, by definition  $\varphi(k)=\prod_{i=1}^m p_i^{a_i-1}(p_i-1)$. As we are assuming that $\ord_k(p)$ is equal to some power of $p$, this implies that $p \mid \varphi(k)$, hence $p \mid (p_i-1)$ for some $i$. But it is clear that this can happens if and only if $p=2$.

Using these considerations, we proved the following result:

\begin{theorem}
Let $f$ be an irreducible polynomial over $\F_p$ of degree $n\ge 2$; then, $f$ satisfies property \eqref{main_property} if and only if either $f(x)=x^{r-1}+ \cdots +1$ with $r$ a prime different from $p$ and $p$ a generator of $(\Z/r\Z)^{\times}$ or $p=2$ and $f$ is a factor of degree $n$ of a cyclotomic polynomial $\Phi_k$ with $n < \varphi(k)$, where $k$ is odd and and $n=\ord_k(2)$ is a power of $2$.
\end{theorem}

\begin{example}
Let us consider the cyclotomic polynomial of order $15$ over $\F_2$; as seen before, we have that  $\Phi_{15}$ factorises into $2$ factors of degree $4$, i.e.
\[
\Phi_{15}(x)=(x^4+x^3+1)(x^4+x+1).
\]
In this case both the factors of $\Phi_{15}$ are polynomials of degree $4$ satisfying property \eqref{main_property}. We point out that if $\Phi_n$ factorises into irreducible factors over $\F_2$ then either none of the factors satisfy the property or all do.
\end{example}

\bibliographystyle{amsalpha}

config: arxiv


TITLE: Prob. 7, Sec. 3.8 in Erwine Kreyszig's INTRODUCTORY FUNCTIONAL ANALYSIS WITH APPLICATIONS: The dual space of a Hilbert space is a Hilbert space.
QUESTION [3 upvotes]: Here's Prob. 7, Sec. 3.8 in Introductory Functional Analysis With Applications by Erwine Kreyszig:

Show that the dual space $H^\prime$ of a Hilbert space $H$ is a Hilbert space with inner product $\langle \ \cdot \ , \ \cdot \ \rangle_1$ defined by
  $$\langle f_z, f_v \rangle_1 \ = \ \overline{\langle z, v \rangle} \ = \ \langle v, z \rangle,$$
  where $f_z(x) = \langle x, z \rangle$ for all $x \in X$.

Now I know that each bounded linear functional $f \in H^\prime$ can be written as $f = f_z$, for a unique $z \in H$ with $\Vert z \Vert = \Vert f \Vert$.
The sum of two bounded linear functionals on any normed space is again a bounded linear functional, and so is the scalar multiple of any bounded linear functional.
Moreover, for each $z \in H$ and for each $w \in H$, we have
$$\left( f_z + f_w \right) (x) = f_z(x) + f_w(x) = \langle x, z \rangle + \langle x, w \rangle = \langle x, z+w \rangle = f_{z+w}(x)$$
for all $x \in H$. So $f_{z+w} = f_z + f_w$.
For $z \in H$ and for any scalar $\alpha$, we have
$$f_{\alpha z} (x) = \langle x, \alpha z \rangle = \overline{\alpha} \langle x, z \rangle = \overline{\alpha} f_z (x) \ \mbox{ for all } \ x \in H.$$
So, $f_{\alpha z} = \overline{\alpha } f_z$.
Thus the mapping $z \mapsto f_z$ of $H$ into $H^\prime$ is surjective, isometric, and conjugate linear.
Moreover, the set of all bounded linear functionals on any normed space is itself a normed space, rather a Banach space.
But how to show that the inner product given by Kreyszig is the (only) natural one (i.e. that which fits into what we already know from the normed space theory)?

REPLY [2 votes]: Let me clarify what the question is asking, since based on your comments you actually already know how to solve it.
As you know, for any complex normed space $(X, \|\cdot\|_X)$, the dual space $X'$ is the complex vector space of all continuous linear functionals $f : X \to \mathbb{C}$, with pointwise addition and scalar multiplication.  We usually equip $X'$ with the "operator norm" $\|\cdot\|_{X'}$ defined by
$$\|f\|_{X'} = \sup_{x \in X, \|x\|_X = 1} |f(x)|.$$
Now in case $X$ is a Hilbert space $(H, \langle \cdot, \cdot \rangle_H)$, we know from the Riesz representation theorem that every $f \in H'$ is of the form $f_z$ for some $z \in H$, where $f_z$ is defined by $f_z(v) = \langle v, z \rangle$.  This suggests constructing an inner product on the complex vector space $H'$: Kreyszig calls this new inner product $\langle \cdot, \cdot \rangle_1$, and it is defined by $\langle f_z, f_v \rangle_1 = \langle v, z \rangle_H$.  You should verify that $\langle \cdot, \cdot \rangle_1$ is an inner product (sesquilinear, Hermitian, positive definite) but this is very easy.
Now, on the complex vector space $H'$ we have a norm $\|\cdot\|_{H'}$ and an inner product $\langle \cdot, \cdot \rangle_1$.  The question is asking you to show that the inner product $\langle \cdot, \cdot \rangle_1$ induces the norm $\|\cdot\|_{H'}$.  That is, you must show that $\|f\|_{H'} = \sqrt{\langle f, f \rangle_1}$ for every $f \in H'$.
If you unwind this, it comes down to showing something like the following: for every $z \in H$, $$\|z\|_H = \sup_{v \in H, \|v\|=1} |\langle v,z \rangle|.$$
Prove this by proving two inequalities.  For one of them, use Cauchy-Schwarz.  For the other, consider an appropriate choice for $v$.
You only have to show this for the specific inner product $\langle \cdot, \cdot \rangle_1$ given in the question, not for an arbitrary inner product on $H'$.  The question is not asking you to show that every inner product on $H'$ induces the norm $\|\cdot\|_{H'}$.  This is trivially false; as a silly example, consider the inner product $\langle \cdot, \cdot\rangle_2$ defined by $\langle f_z, f_v \rangle_2 = 2 \langle v, z \rangle_H$.  If you use the axiom of choice, you can produce other inner products on $H'$ that induce completely different (and uninteresting) topologies.