# Integrated Decision Gradients: Compute Your Attributions Where the Model
Makes Its Decision

Chase Walker
Department of ECE
University of Central Florida
Orlando, FL, USA
<EMAIL_ADDRESS>
&Sumit Jha
Computer Science Department
University of Texas at San Antonio
San Antonio, TX, USA
<EMAIL_ADDRESS>
&Kenny Chen
Lockheed Martin Corporation
Orlando, FL, USA
<EMAIL_ADDRESS>
&Rickard Ewetz
Department of ECE
University of Central Florida
Orlando, FL, USA
<EMAIL_ADDRESS>

###### Abstract

Attribution algorithms are frequently employed to explain the decisions of
neural network models. Integrated Gradients (IG) is an influential attribution
method due to its strong axiomatic foundation. The algorithm is based on
integrating the gradients along a path from a reference image to the input
image. Unfortunately, it can be observed that gradients computed from regions
where the output logit changes minimally along the path provide poor
explanations for the model decision, which is called the _saturation effect_
problem. In this paper, we propose an attribution algorithm called integrated
decision gradients (IDG). The algorithm focuses on integrating gradients from
the region of the path where the model makes its decision, i.e., the portion
of the path where the output logit rapidly transitions from zero to its final
value. This is practically realized by scaling each gradient by the derivative
of the output logit with respect to the path. The algorithm thereby provides a
principled solution to the saturation problem. Additionally, we minimize the
errors within the Riemann sum approximation of the path integral by utilizing
non-uniform subdivisions determined by adaptive sampling. In the evaluation on
ImageNet, it is demonstrated that IDG outperforms IG, left-IG, guided IG, and
adversarial gradient integration both qualitatively and quantitatively using
standard insertion and deletion metrics across three common models.

##  1 Introduction

The access to internet-scale data and compute power has fueled the success of
black box neural network models for applications such as disease detection
[1], image synthesis [2], and protein folding [3]. The phenomenal performance
of these networks comes from the large number of parameters and non-linear
interactions among them. The complex and high dimensional dynamics makes it
difficult to understand and visualize why a neural network makes a particular
decision. To establish trustworthiness in neural network models, noteworthy
research efforts have been devoted to interpretability and explainability [4].
Attribution methods provide model explanation by computing the contribution of
each input feature to a model decision. Attribution methods broadly fall into
perturbation based methods [5, 6], backpropagation based methods [7, 8], and
gradient based methods [9, 10]. Gradient based methods are promising due to
their strong axiomatic foundation, and model-agnostic implementation [10].

Gradient based methods compute attribution maps by capturing the gradients at
the model inputs with respect to the model outputs [9]. However, gradients
computed with respect to important input pixels may be zero due to the non-
linear activation functions. Integrated Gradients (IG) solved this problem by
integrating the gradients along a path from a baseline reference image to the
input image. Unfortunately, it can be observed that gradients from regions of
the path where the output logit changes minimally (e.g. is saturated) provide
poor explanations for the model decision [11]. This phenomena is called the
_saturation effect_ problem. Solution templates to solve the saturation
problem include: selecting non-straight-line paths [12, 13], path truncation
[11], and post processing methods that use thresholding [14] as well as
averaging across blurred inputs [15]. While these methods improve attribution
quality, they do not provide a principled solution to the saturation problem.

Figure 1: (left) An overview of the adaptive sampling algorithm, and the IDG
attribution method. (right) A preliminary visual comparison of IDG with IG
[10], LIG [11], GIG [12], and AGI [13].

In this paper, we propose a new attribution method called Integrated Decision
Gradients (IDG). We call the portion of the path where the output logit
rapidly transitions from zero to its final value the _decision region_. The
IDG algorithm focuses on integrating gradients from the decision region of the
path integral. This is practically realized by scaling each gradient by the
derivative of the output logit with respect to the path. The scaling factor
rewards gradients in the decision region and penalizes gradients from outside
the decision region. The main contributions of this paper are summarized as
follows:

  * •

We propose a new attribution method called IDG that satisfies a path integral
sensitivity axiom and provides a principled solution to the saturation
problem.

  * •

We present an adaptive sampling technique to select non-uniform subdivisions
for the Riemann approximation of the path integral. The non-uniform
subdivisions reduce computational errors (or runtime overheads) compared with
using uniform subdivisions.

  * •

Compared with IG [10], Left-IG (LIG) [11], Guided IG (GIG) [12], and
Adversarial Gradient Integration (AGI) [13], IDG improves both the qualitative
and quantitative results.

The remainder of the paper is organized as follows: related work is examined
in Section 2, the IDG attribution method in Section 3, the adaptive sampling
algorithm is proposed in Section 4, experimental evaluation is presented in
Section 5, and the paper is concluded in Section 6.

##  2 Related Work

In this section, we first review the limitations of directly using gradients
as attributions. Next, we review integrated gradients and assess the
saturation effect problem within path integrals.

###  2.1 Limitations of Using Gradients as Attributions

Attributions are defined to be the contribution of each input feature to the
model output decision. An attribution method satisfies the axiom of
_sensitivity_ if a single feature that differs between a baseline and input -
which produce different output predictions - is given a non-zero attribution.
Additionally, if a neural network is not affected by changing a variable, then
that variable’s attribution shall be zero [10]. Computing the gradient of the
inputs with respect to the output logit is a promising method for computing
attributions [9]. However, the use of non-linear activation functions causes
the sensitivity axiom to be violated [10], which is shown in Example 1 below.

######  Example 1.

Consider a function $F=1-ReLU(1-x)$, a baseline $x^{\prime}=0$, and an input
$x=2$. For $x^{\prime}=0$, the function $F$ is equal to $0$, and for $x=2$,
the function $F$ is equal to $1$. Since changing $x$ from $0$ to $2$ affects
the output of $F$, the attribution w.r.t. the feature $x$ should be non-zero.
However, $\partial F/\partial x=0$ at $x=2$, which results in an attribution
of $0$ [10].

Integrated gradients offers a solution to computing attributions that
satisfies the sensitivity axiom.

Figure 2: The figure illustrates the IG attribution method and saturation
effects within path integrals. The top row shows interpolated inputs and the
second row shows the corresponding gradients. The IG attribution map (shown to
the right) is the average of the gradients. The third row shows the
logit-$\alpha$ curve, which defines the decision and saturation regions. It
can be observed that the gradients from the decision region are of higher
quality than the saturation region.

###  2.2 Integrated Gradients

Integrated Gradients computes attributions by integrating gradients on a
straight line between a reference image and an input image [10]. Let $F$ be
the function realizing the output logit of interest. $IG_{i}$ with input image
$x$ is mathematically defined using a path-integral [10], as follows:

$IG_{i}(x)=(x_{i}-x^{\prime}_{i})\times\int_{\alpha=0}^{1}\frac{\partial
F(x^{\prime}_{i}+\alpha\times(x_{i}-x^{\prime}_{i}))}{\partial x_{i}}d\alpha$
(1)

where $x^{\prime}$ is a black baseline image, $\alpha\in[0,1]$ parameterizes
the straight-line path between $x^{\prime}$ and $x$, $x_{i}$ and
$x^{\prime}_{i}$ represent a single pixel of their respective images, and
$IG_{i}$ is therefore the attribution of pixel $i$ in the input image.

The IG attribution method is illustrated in Figure 2. The top row shows
interpolated inputs, the second row shows the corresponding input gradients,
the third row visualizes the output logit with respect to the path. The IG
attribution map is equal to the sum of the gradients in the second row. The
use of a path-integral ensures that gradients from regions of $F$ where
$\partial F/\partial x_{i}$ is non-zero are computed. In Example 1, IG will
compute gradients from the region $[0,$ $1]$, where $\partial F/\partial x=1$.
The resulting attribution w.r.t. $x$ is $2$, i.e., the attribution is non-zero
and sensitivity is satisfied. Nevertheless, many attributions computed using
IG are still noisy due to saturation effects [11].

###  2.3 Saturation Effects within Path-Integrals

To introduce and understand the _saturation effect_ problem within path-
integrals, we examine the performance of the IG attribution method in in
Figure 2. We study the quality of the computed gradients with respect to the
decision and saturated regions of the path integral. It can be observed that
(i) gradients from the saturation regions are of low quality and (ii)
gradients from the decision region are of high quality. The conclusion is
rather straight forward to understand. If the model output does not increase
while moving $\triangle\alpha$ along the path, it is intuitive that the
corresponding gradients are not important to the model decision. Conversely,
if the output logit changes rapidly while moving $\triangle\alpha$ along the
path, those gradients have a strong impact on the model decision.

This raises the rudimentary question: Is it possible to design a path integral
that focuses on computing gradients from the region where the model decision
is made and the highly informative gradients are located? It can for example
be observed in Figure 2 that the gradients computed at $\alpha=0.02$ alone
provide an excellent explanation for the model decision.

##  3 Integrated Decision Gradients

In this section, we propose a new attribution method called Integrated
Decision Gradients (IDG). We outline the motivation behind the design of IDG,
explain the concept of importance factors, and provide the definition as well
as a visualization of IDG.

###  3.1 Motivation

Path integrals integrate gradients from a reference image to an input target
image. A fundamental challenge is to determine the ideal importance of each
gradient. Based on the analysis in the previous section, we define a new
_sensitivity axiom_ for path integrals. Next, we introduce the concept of an
importance factor, which is used to construct an attribution algorithm that
satisfies the axiom.

#### Axiom: Sensitivity (path integrals)

Let $F$ be the output of a neural network. For every point within a path
integral parameterized by a parameter $\alpha$, an attribution method
satisfies Sensitivity (path integrals) if there is no contribution to the
attribution result when $\partial F/\partial\alpha$ is equal to zero. If
$\partial F/\partial\alpha$ is non-zero, there should be a non-zero
contribution to the attribution result.

Figure 3: This figure illustrates the relationship between importance factor
magnitude and gradient quality. Higher importance factors are directly related
to higher quality gradients.

None of the existing attribution methods based on path integrals satisfy the
axiom [10, 11, 12, 13]. The traditional IG method places an equal weight on
all gradients [10], even those that occur in the saturation region where
$\partial F/\partial\alpha=0$. The Left-IG attribution attempts to solve this
by truncating the path integral after the output logit has reached $90\%$ of
its final value [11]. This assigns a weight of zero and one to gradients from
the approximate saturation and decision regions respectively, which does not
guarantee that the axiom is satisfied. GIG and AGI use non-straight line paths
that attempt to avoid integrating gradients from saturated regions [12, 13],
which does also not guarantee that the Sensitivity (path integrals) axiom is
satisfied.

To satisfy the axiom, we conjecture that the importance of each gradient
should be proportional to the impact on the model output, which is
conceptually shown in Figure 3. Inspired by this, we define an _importance
factor_ , as follows:

$\displaystyle IF(\alpha)=\frac{\partial
F(x^{\prime}+\alpha(x-x^{\prime}))}{\partial\alpha}$ (2)

where $IF(\alpha)$ is the importance of the gradient computed at $\alpha$.
Next, we define an attribution method that satisfies the Sensitivity (path
integrals) axiom based on scaling each gradient with the importance factor in
Eq (2).

###  3.2 Definition of Integrated Decision Gradients

In this subsection, we formally define the IDG attribution algorithm. Given a
neural network represented by function $F:R^{n}\to[0,1]$, an input vector $x$,
and given $F$ exists over the range $\alpha\in[0,1]$, IDG assigns an
importance factor to each input feature $x_{i}$ with respect to the model
output, using the following equation:

$IDG_{i}(x)=\underbrace{(x_{i}-x^{\prime}_{i})\times\int_{\alpha=0}^{1}\frac{\partial
F(x^{\prime}_{i}+\alpha(x_{i}-x^{\prime}_{i}))}{\partial
x_{i}}}_{\text{Traditional IG}}\times\underbrace{\frac{\partial
F(x^{\prime}_{i}+\alpha(x_{i}-x^{\prime}_{i}))}{\partial\alpha}}_{\text{Importance
Factor}}d\alpha$ (3)

The IDG attribution method is equivalent to IG in Eq (1) but with each
gradient scaled with the importance factor in Eq (2). The importance factor is
equivalent to the derivative of the logit-$\alpha$ curve in the bottom of
Figure 2. The importance factors scale-up gradients from the decision region
and scale-down gradients from saturated regions, respectively. Therefore, IDG
provides a principled solution to the saturation problem, and satisfies the
Sensitivity (path integrals) axiom by definition.

The path integral is practically computed using the Riemann sum approximation
[10], as follows:

$IDG_{i}(x)=(x_{i}-x^{\prime}_{i})\times\frac{1}{m}\times\sum_{k=1}^{m}\frac{\partial
F(x^{\prime}_{i}+\frac{k}{m}\ \times(x_{i}-x^{\prime}_{i}))}{\partial
x_{i}}\times\frac{\partial F(x^{\prime}_{i}+\frac{k}{m}\
\times(x_{i}-x^{\prime}_{i}))}{\partial\alpha}$ (4)

where $m$ is the number of steps for approximation. We will further discuss
the selection of the step size and its impact on the approximation error in
Section 4.

We illustrate IDG with an example in Figure 4. First, looking at the left side
of the figure, the top row shows the logit-$\alpha$ curve associated with the
input image. The second row shows the derivative of this curve, i.e.,
$\partial F/\partial\alpha$ in Eq (2). The third row shows the interpolated
inputs for selected alpha values and the fourth row shows the gradients
computed by IG for these inputs. The last row visualizes the effect of IDG by
scaling the gradients above by the importance factors from the second graph.
The importance factors scale up the magnitude of the gradients from the
decision region while scaling down the magnitude of the gradients from the
saturated regions. In the figure, it can be observed that, in particular, the
attributions from $\alpha=0.005$ are scaled up. On the right of the figure, we
show the original image, and the attributions generated by IG, LIG, GIG, AGI,
and IDG. The attributions computed using IDG are substantially less noisy than
all competitors. We note that GIG has a low amount of noise, but IDG has more
focused attributions on the highlighted features.

Figure 4: A full visualization of how IDG uses importance factors to eliminate
saturation effects. The top row shows the logit-$\alpha$ curve. The next row
shows the derivative of the curve, i.e., the importance factors with respect
to $\alpha$. The third row shows the interpolated images, the fourth shows the
associated gradients, and the bottom row shows these gradients scaled with the
corresponding importance factors. The right side shows the input image, and
the attributions computed using IG [10], LIG [11], GIG [12], AGI [13], and
IDG.

##  4 Adaptive Sampling Algorithm

In this section, we first analyze the errors within the Riemann sum
approximation of the IDG path integral for uniform subdivisions. Next, we
propose an adaptive sampling technique to minimize the approximation errors
using non-uniform subdivisions. In the supplementary results, we show that the
adaptive sampling only results in major improvements for IDG. The impact of
the adaptive sampling on regular IG is minor.

Figure 5: This figure shows the motivation for the adaptive sampling
algorithm. The image (a) is the input to the attributions in the figure. The
graph (b) demonstrates how the attribution error decreases as step count
increases. Columns (c), (d), and (e) of attributions and graphs show the
relationship between sample locations and IDG quality as $50$, $250$, and
$600$ steps are used respectively. We show that as the number of steps
increases, the quality of IDG grows greatly, influencing the adaptive sampling
algorithm. Lastly, column (f) shows the equivalent result of column (e)
achieved by using adaptive sampling with 50 steps.

###  4.1 Motivation

The errors within the Riemann approximation of the IDG path integral can be
calculated, as follows:

$\epsilon(x,n)=\lim_{m\to\infty}IDG_{i}(x,m)-IDG_{i}(x,n)$ (5)

where $\epsilon(x_{i},n)$ is the approximation error for attribution $x_{i}$
when computing the integral with $n$ uniform subdivisions. $n$ and $m$ are the
number of steps used within the Riemann sum approximation in Eq (4).

We analyze the approximation error and the impact on the attributions in
Figure 5. The graph (b) shows the average error across all the pixels in the
attribution map with respect to the number of used steps $n$. Since a low step
count results in a lack of samples in the decision region, a large number of
steps are required for a good approximation. The image (a) is the input for
the four columns (c), (d), (e), and (f) of attributions and graphs. The
columns show the quality of the attributions with respect to the number of
steps and type of subdivision. It is observed from the graphs that taking more
samples in the decision region greatly improves IDG attribution quality.
Therefore, to obtain high IDG quality without a prohibitive number of steps,
we design a new adaptive sampling algorithm - seen in Figure 5 (f) - that uses
non-uniform subdivisions concentrated on the decision region.

###  4.2 Adaptive Sampling Methodology

It is desirable to sample the high quality gradients that lie in the decision
region to improve the quality of the attained attributions. In Algorithm 1, we
show how the adaptive sampling algorithm is used with IDG. Our approach is
based on first pre-characterizing the logit-$\alpha$ curve with $N$ uniform
subdivisions in lines 3 - 7. Next, $M$ subdivisions are non-uniformly
distributed within the $N$ regions based on logit growth and IDG is calculated
in lines 8 - 15. Because there are $M$ total samples, line 11 executes
$O(N+M)$ times. In practice it is best if $N=M$ (this is shown in the
supplementary materials) therefore the algorithm runtime is $O(N)$.

As seen in Figure 5 (e) and (f), combining this adaptive sampling algorithm
with IDG creates attributions as strong as IDG with 600 steps while only using
50 steps. Figure 1 provides a high-level overview of this new IDG process. The
figure shows that when given an input image and a number of steps, the
adaptive sampling algorithm calculates non-uniform subdivisions based on logit
growth. These are then used as input for IDG where the gradient at each
location is calculated and then weighted, producing the final attribution. In
this figure, the IDG sampling graph shows that $31$ out of $50$ samples are
placed in the decision region $\alpha\in[0.0,0.2]$, where the logit changes
from $0$ to $7.2$.

Input: Model $F$, image $x$, baseline $x^{\prime}$, pre-characterization steps
$N$, number of IDG steps $M$
Output: An attribution map $A$

Algorithm 1 Computing IDG with Adaptive Sampling

1: $x^{0}=x^{\prime}$

2: $x^{N-1}=x$

3: $samples[0]=0$ // Pre-characterization of logit-$\alpha$ curve

4: for $i=0$ to $N-1$ do

5: $x^{i+1}=x^{\prime}+\frac{i}{N}\times(x-x^{\prime})$

6: $samples[i+1]=round(\frac{F(x^{i+1})-F(x^{i})}{F[x^{N-1}]-F[x^{0}]}\times
M$)

7: end for// Computation of IDG with non-uniform samples

8: for $i=0$ to $N$ do

9: for $j=0$ to $samples[i]$ do

10: $\alpha=\frac{i}{N}+\frac{j}{N\times samples[i]}$

11: $x^{i}=x^{\prime}+\alpha\times(x-x^{\prime})$

12: $IDG[i]=\frac{\partial F(x^{i})}{\partial x}\times\frac{\partial
F(x^{i})}{\partial\alpha}\times\frac{1}{N\times samples[i]}$

13: end for

14: end for

15: return $A=mean(IDG)$

##  5 Experimental Results

In this section, we will evaluate the effectiveness of the proposed method. We
perform our experiments in PyTorch using the 2012 validation set of ImageNet
[16] on NVIDIA A40 GPUs. According to ML CO2 impact, the experimental
evaluation released $43.6$ kg of CO2 with zero offset [17]. The attributions
computed using Algorithm 1 are called IDG. We compare our method with IG [10],
left-IG [11], guided IG [12], and adversarial gradient integration [13]. We
use Captum for the implementation of IG, whereas left-IG, GIG, and AGI are
taken from their respective repositories [18, 19, 20, 21]. We evaluate the
quality of the computed attributions both quantitatively and qualitatively.

In Table 1, we quantitatively evaluate the attributions using standard
perturbation testing which measures the importance of the pixels in an
attribution via an area under the curve (AUC) score. A total of four tests are
presented with three insertion methods and one deletion method from the
authors of RISE and XRAI [22, 14] which are described in Section 5.1. The
table compares the computed attribution quality for the first $5000$ images of
the ImageNet dataset such that five images are taken from each of the $1000$
classes. The five attribution methods are evaluated with three models trained
on ImageNet. We selected ResNet101 (R101), ResNet152 (R152), and ResNeXt
(RNXT) as pre-trained models from PyTorch and use the newest ImageNet weights
available (V2 for the ResNet models and V1 for ResNeXt) [23, 24, 25].

Qualitatively, we present a subset of five examples in Figure 6 which are
gathered gathered using the ResNet101 model and the method parameters
explained below. We provide a larger selection of examples in the
supplementary materials for visual comparison.

Inputs are reshaped to ($224$, $244$) for all three presented models. This
image processing follows the attribution documentation provided by Captum
[18]. The RISE, AIC, and SIC tests use the default parameters found from their
respective repositories [26, 27]. The IG and LIG attribution methods use $50$
steps and a black baseline image. GIG uses the default parameters found at
[20]. AGI uses the default parameters found at [21]. Lastly, IDG is used with
50 steps and a black baseline image. For all the methods, we use a single
baseline only.

###  5.1 Quantitative Evaluation Metrics

The evaluation metrics are built upon the intuition that the highest
attribution values should correspond to those features that contribute more to
the classification of the target class [22, 14]. The process starts from the
most important pixels and starts deleting (inserting) them from the original
image (to a blurred image for insertion) until only a black (the original)
image remains. At each step, the softmax score (or accuracy) is calculated.
This gives us an ROC curve from base image to final image, which is used to
compute the AUC score for a given attribution. This AUC value is computed for
each image and then averaged out over the entire test data selection. For the
insertion game, a higher AUC score indicates a better attribution and for the
deletion game, a lower AUC score indicates better performance. The two sets of
methods presented from Petsiuk, et al and Kapishnikov, et al. take different
approaches to the insertion process [22, 14].

In RISE, the insertion (deletion) test which starts (ends) with a Gaussian
blurred (black) image [22]. In their implementation pixels are added (deleted)
in equal amounts during the test process. Given an NxN image, the test will
change the image by N pixels at a time over N steps.

Kapishnikov, et al. present the Accuracy Information Curve (AIC) and Softmax
Information Curve (SIC) in their XRAI paper [14]. The AIC test gives each
perturbation step a score of $0$ or $1$ for an incorrect or correct
classification and SIC uses softmax as previously discussed. For pixel
perturbation, these methods use a schedule that non-linearly removes groups of
pixels from the image in increasingly large amounts. The last difference from
the RISE insertion test is the blurring method, where the initial image is now
blurred in segments, each having its own noise distribution.

Table 1: Comparison of attributions using the AIC, SIC, insertion, and deletion tests Metric | Model | IG [10] | LIG [11] | GIG [12] | AGI [13] | IDG
---|---|---|---|---|---|---
| R101 | 0.571 | 0.589 | 0.626 | 0.675 | 0.701
AIC ($\uparrow$) | R152 | 0.575 | 0.616 | 0.646 | 0.686 | 0.718
| RNXT | 0.580 | 0.611 | 0.634 | 0.654 | 0.730
| R101 | 0.498 | 0.522 | 0.559 | 0.609 | 0.638
SIC ($\uparrow$) | R152 | 0.508 | 0.552 | 0.582 | 0.619 | 0.659
| RNXT | 0.478 | 0.518 | 0.532 | 0.554 | 0.620
| R101 | 0.498 | 0.535 | 0.547 | 0.561 | 0.592
Insertion ($\uparrow$) | R152 | 0.517 | 0.562 | 0.565 | 0.577 | 0.615
| RNXT | 0.276 | 0.299 | 0.296 | 0.307 | 0.324
| R101 | 0.181 | 0.148 | 0.155 | 0.172 | 0.108
Deletion ($\downarrow$) | R152 | 0.202 | 0.148 | 0.164 | 0.190 | 0.118
| RNXT | 0.101 | 0.078 | 0.082 | 0.104 | 0.068

###  5.2 Comparison With Previous Work

In Table 1, attribution quality is evaluated using the AIC and SIC insertion
metrics and the RISE insertion and deletion metrics. We use an arrow to denote
if larger (arrow up) or smaller (arrow down) scores are better. The best score
for each model and test type is in bold. Additionally we provide how many
times a given method outperforms all other methods in the last row of the
table.

It can be observed in Table 1 that IDG achieves a consistent improvement over
IG, LIG, GIG, and AGI across all twelve of the tests presented. Comparing IDG
to IG and LIG clearly indicates the ability of IDG to mitigate saturation
effects in path-based methods while retaining the most important gradient
information. When compared to AGI and GIG, the large margin of improvement in
the scores shows that IDG presents a more complete solution to the saturation
problem than these methods. Overall, IDG outperforms all of the attribution
methods in the comparison, achieving new state-of-the-art performance.

For qualitative analysis, we compare IG, LIG, GIG, AGI, and IDG in Figure 6.
All attributions are computed as previously described. The comparison is
performed using images of a “Guenon", “Submarine", “Tripod", “African Hunting
Dog", and “Warplane" taken from ImageNet [16].

Across the five selections, IDG clearly produces much sharper attributions
than IG and LIG, further verifying that it solves the saturation problem
present in these methods. When compared to GIG, IDG also has superior
performance in all of the images. For the Tripod example, even though GIG has
relatively low noise, IDG has stronger attributions on the tripod in the
foreground, and the one in the background as well. Lastly, when comparing to
AGI, it can be seen AGI generally has low extraneous noise in the
attributions. However, IDG provides tighter, and sharper attributions on the
class subject in the images, therefore the results are better.

The images clearly show that IDG improves visual quality over the other
methods. IDG generates attributions with less random noise, and sharper
attributions, showing its ability to solve the saturation problem.
Additionally it shows its ability to outperform the methods which use non-
straight-line paths. We provide an additional 50 visual comparisons in the
supplementary results section.

Figure 6: Qualitative comparison of attributions computed using the IG [10],
LIG [11], GIG [12], and AGI [13], and IDG methods. It is seen that IDG solves
the saturation problem and outperforms the state-of-the-art path-based
attribution methods in visual quality.

##  6 Discussion

In this paper, we propose a new attribution method called Integrated Decision
Gradients (IDG). The key idea of IDG is to perform the path integral while
weighting sampled gradients by their associated logit growth. This amplifies
gradients located in the decision region, and negates those from the
saturation region, solving the saturation issue. In contrast, traditional IG
integrates gradients between the same images while giving all gradients equal
weight, saturated or not, causing the majority of saturated gradients to
dominate the output. Additionally, we provide evidence that the decision
region of the path integral is where the best gradients lie. With this, we
present an adaptive sampling algorithm which densely samples the decision
region without runtime penalty, improving IDG performance. We show
qualitatively and quantitatively that IDG reaches state-of-the-art performance
in the path-based attribution field. In our future work, we plan to apply IDG
concepts to other attribution methods to further enhance attribution quality.
We also plan to employ IDG within practical real-world applications. The code
to replicate the results presented in this paper is available at:
https://github.com/chasewalker26/Integrated-Decision-Gradients.

#### Limitations

We present quantitative and qualitative results that show the proposed method
outperforms those in its field. However, there does not yet exist criteria to
perfectly examine what makes a good attribution. Therefore, we provide the
best, currently accepted evaluation of our method.

## Acknowledgements

This work was partly supported by the Lockheed Martin University Engagement
Program, the Florida High Tech Corridor Matching Grants Program, and the DARPA
cooperative agreement #HR00112020002. The views, opinions and/or findings
expressed are those of the authors and should not be interpreted as
representing the official views or policies of the Department of Defense or
the U.S. Government.

## References

  * [1] Meherwar Fatima, Maruf Pasha, et al.  Survey of machine learning algorithms for disease diagnostic.  Journal of Intelligent Learning Systems and Applications, 9(01):1, 2017.
  * [2] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695, 2022.
  * [3] Milot Mirdita, Konstantin Schütze, Yoshitaka Moriwaki, Lim Heo, Sergey Ovchinnikov, and Martin Steinegger.  Colabfold: making protein folding accessible to all.  Nature Methods, pages 1–4, 2022.
  * [4] Arun Das and Paul Rad.  Opportunities and challenges in explainable artificial intelligence (xai): A survey.  arXiv preprint arXiv:2006.11371, 2020.
  * [5] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In European conference on computer vision, pages 818–833. Springer, 2014.
  * [6] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  "why should i trust you?" explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135–1144, 2016.
  * [7] Jost Tobias Springenberg, Alexey Dosovitskiy, Thomas Brox, and Martin Riedmiller.  Striving for simplicity: The all convolutional net.  arXiv preprint arXiv:1412.6806, 2014.
  * [8] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In Proceedings of the IEEE international conference on computer vision, pages 618–626, 2017.
  * [9] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps, 2013.
  * [10] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 3319–3328. JMLR.org, 2017.
  * [11] Vivek Miglani, Narine Kokhlikyan, Bilal Alsallakh, Miguel Martin, and Orion Reblitz-Richardson.  Investigating saturation effects in integrated gradients, 2020.
  * [12] A. Kapishnikov, S. Venugopalan, B. Avci, B. Wedin, M. Terry, and T. Bolukbasi.  Guided integrated gradients: an adaptive path method for removing noise.  In 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5048–5056, Los Alamitos, CA, USA, jun 2021. IEEE Computer Society.
  * [13] Deng Pan, Xin Li, and Dongxiao Zhu.  Explaining deep neural network models with adversarial gradient integration.  In Zhi-Hua Zhou, editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, pages 2876–2883. International Joint Conferences on Artificial Intelligence Organization, 8 2021.  Main Track.
  * [14] A. Kapishnikov, T. Bolukbasi, F. Viegas, and M. Terry.  Xrai: Better attributions through regions.  In 2019 IEEE/CVF International Conference on Computer Vision (ICCV), pages 4947–4956, Los Alamitos, CA, USA, nov 2019. IEEE Computer Society.
  * [15] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise, 2017.
  * [16] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei.  ImageNet Large Scale Visual Recognition Challenge.  International Journal of Computer Vision (IJCV), 115(3):211–252, 2015.
  * [17] Alexandre Lacoste, Alexandra Luccioni, Victor Schmidt, and Thomas Dandres.  Quantifying the carbon emissions of machine learning.  arXiv preprint arXiv:1910.09700, 2019.
  * [18] Narine Kokhlikyan, Vivek Miglani, Miguel Martin, Edward Wang, Bilal Alsallakh, Jonathan Reynolds, Alexander Melnikov, Natalia Kliushkina, Carlos Araya, Siqi Yan, et al.  Captum: A unified and generic model interpretability library for pytorch.  arXiv preprint arXiv:2009.07896, 2020.
  * [19] Vivek Miglani, Narine Kokhlikyan, Bilal Alsallakh, Miguel Martin, and Orion LLC Reblitz-Richardson.  Left-IG Code Repository, 2020.  Available at https://github.com/vivekmig/captum-1/tree/ExpandedIG.
  * [20] A. Kapishnikov, S. Venugopalan, B. Avci, B. Wedin, M. Terry, and T. Bolukbasi.  Gig code repository, 2021.  Available at https://github.com/PAIR-code/saliency/tree/master/saliency/core.
  * [21] Deng Pan, Xin Li, and Dongxiao Zhu.  Agi code repository, 2021.  Available at https://github.com/pd90506/AGI.
  * [22] Vitali Petsiuk, Abir Das, and Kate Saenko.  Rise: Randomized input sampling for explanation of black-box models.  In Proceedings of the British Machine Vision Conference (BMVC), 2018\.
  * [23] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [24] Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He.  Aggregated residual transformations for deep neural networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1492–1500, 2017.
  * [25] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Proceedings of the 33rd International Conference on Neural Information Processing Systems. Curran Associates Inc., 2019.
  * [26] Vitali Petsiuk, Abir Das, and Kate Saenko.  Rise code repository, 2018.  Available at https://github.com/eclique/RISE.
  * [27] A. Kapishnikov, S. Venugopalan, B. Avci, B. Wedin, M. Terry, and T. Bolukbasi.  Aic code repository, 2019.  Available at https://github.com/PAIR-code/saliency/tree/master/saliency/metrics.
  * [28] Content European Commission, Directorate-General for Communications Networks and Technology.  Ethics guidelines for trustworthy ai, 2019.
  * [29] Mark Ryan.  In ai we trust: ethics, artificial intelligence, and reliability.  Science and Engineering Ethics, 26(5):2749–2767, 2020.

##  Appendix A Appendix

In this appendix we provide additional information that did not fit in the
bounds of the paper. In Section A.1 we provide detailed analysis of the
selection of $N$ and $M$ for the adaptive sampling algorithm. In Section A.2
we provide further explanation of the impact of the AS algorithm, showing that
IDG provides the true solution to the saturation problem. In Sections A.3 and
A.4 we provide information on the licensing of libraries used in experimental
evaluation and we discuss the potential ethical impact of the proposed method.
Lastly, in Section A.5 we present 50 additional qualitative visual comparisons
of our proposed method against those presented in the manuscript.

###  A.1 Ablation Study for Adaptive Sampling

The adaptive sampling algorithm has two parameters $N$ and $M$. $N$ is the
number of samples used in the pre-characterization of the logit-$\alpha$
curve. $M$ is the number of samples used in the computation of IDG using non-
uniform subdivisions. Three types of selections of $N$ and $M$ are possible:
$N<M$, $N==M$, and $N>M$. Assuming $M$ or $N$ is set to $50$, which is a
common step count for path-based methods, we provide analysis of which
selection provides the best result via an ablation study.

In Figure 7, we present an ablation study on the selection of $N$ and $M$. In
(a), $M$ is set to $50$ and we take the average of the deletion score [22]
over $10$ images as $N$ is varied from $5$ to $100$ by increments of $5$. In
(b), $N$ is set to $50$ and the deletion scores are gathered as before where
$M$ is varied instead. We see from graph (a) that low values of $N$ produce
poor results and the transition from $5$ to $20$ results in a large drop in
deletion score. We see a similar case in (b) where the score improves as $M$
increases. We note that stable performance is seen on both graphs where
$N==M$.

Figure 7: The change in the deletion score of IDG with AS averaged over $10$
images by varying (a) $N$ and (b) $M$. In the graphs, $N$ ($M$) is varied from
$5$ to $100$ while $M$ ($N$) is set to $50$. It is seen that the most stable
scores are located where $N==M$.

We conclude that selecting $N==M$ results in proper estimation of the
importance factors given the $M$ IDG steps available for placement. While
selecting $N<M$ may provide equally strong results, it may provide poor
results without meaningful runtime improvement, therefore $N==M$ is chosen.
Additionally, we note selecting an $N>M$ does not improve the score enough for
the associated runtime penalty.

###  A.2 IDG Is the Solution to the Saturation Problem

In the manuscript we present adaptive sampling as a method to improve the
ability of the proposed IDG method to solve the saturation problem. Adaptive
sampling takes advantage of the importance factors to perform non-uniform
sampling that focuses on the region of growth. This provides a large ratio of
high quality gradients to saturated gradients from which IDG can generate an
attribution. However, AS alone is not a solution to the saturation problem,
which we will demonstrate by evaluating IG with adaptive sampling.

In Figure 8, for the given input image we compare the attributions generated
by IG with uniform sampling (US), IG with adaptive sampling, IDG with uniform
sampling, and IDG with adaptive sampling. When comparing IG with US and IG
with AS, we see a small reduction in noise in the AS attribution, as there are
inherently less saturated gradients captured when AS is applied to IG.
However, due to IG equally weighting all gradients, the saturated gradients
still dominate the output, illustrating that AS alone cannot solve the
saturation problem. However, when viewing IDG with US compared to both IG
attributions, we see a vast improvement to attribution quality, illustrating
IDG’s ability to solve the saturation problem. Furthermore, when AS is is
applied to IDG, its ability to solve the saturation problem increases.

Figure 8: Given the input, attributions created by IG using uniform sampling
(US) and AS are compared to attributions created by IDG with US and AS. Since
AS applied to IG does not meaningfully improve performance over IG with US,
and IDG with US provides a higher quality attribution, we determine that IDG,
not AS, is the solution to the saturation problem. This is further exemplified
by the improvement seen in IDG with AS, reinforcing the idea that AS gives IDG
access to better gradients, but IDG is the solution to the saturation problem.

As adaptive sampling does not meaningfully improve IG performance and IDG with
US provides much stronger attributions than IG, we verify that integrated
decision gradients is the solution to the saturation problem. We reiterate
that adaptive sampling is used to provide IDG access to gradients of a higher
quality than US does, therefore improving its performance, but not acting as
the solution to the saturation problem.

###  A.3 Licenses of Use

The 2012 validation set of ImageNet [16] is used under the BSD 3-Clause
License. The insertion, deletion, AIC, and SIC are used under the MIT and
Apache 2.0 licenses respectively [22, 14]. The IG [10] and LIG [11]
attribution methods as well as the PyTorch repository [25] are used with the
BSD 3-Clause license, GIG [12] is used under Apache License 2.0, and AGI [13]
is used under the MIT license.

###  A.4 Broader Impact

While explainable AI endeavours to increase human trust in AI systems, there
is debate about what trust of AI _is_ , and if unjustified trust can be
harmful. A prominent proponent of complete AI trust is the European
Commission’s High-level Expert Group on AI (HLEG) [28]. This promotion of this
trust is seen as harmful by some who believe creating a reliance on AI will
have a negative impact on humanity and that AI can never truly be trusted due
to its nature [29]. Explainable AI, as presented in this paper, is intended to
allow better validation and understanding of models which are in use or
proposed for use.

###  A.5 Additional Visual Comparisons

To validate the quantitative performance presented in the paper, we visually
compare IG [10], LIG [11], GIG [12], AGI [13], and IDG with a larger number of
examples. We present 50 example images on pages 14 - 18. There are six columns
per example. From left to right the columns are: the input image, IG, LIG,
GIG, AGI, and IDG. These labels are provided above the columns on each page
and the class of the input image is provided to its left. The images are from
the ImageNet validation set. We use ResNet101 on pages 14 - 16, ResNet152 on
page 17, and ResNeXt on page 18 [23, 24]. The attributions are generated with
the same parameters as the quantitative testing.

The presented attributions are analyzed visually. A stronger attribution is
defined by reduction of noise in areas irrelevant to the object of the image,
and stronger attribution (darker color) in areas where the object exists.
After visual analysis, we believe IDG presents a sharper attribution than all
of the methods presented for a majority of the provided examples. This
thorough qualitative analysis provides further proof of the strength of the
proposed IDG method.

See pages 1-3 of IDG_supplementary.pdf See pages - of
IDG_supplementary_R152.pdf See pages - of IDG_supplementary_RNXT.pdf


#  Melnikov Method for Perturbed Completely Integrable Systems 00footnotetext:
2020 Mathematics Subject Classification. Primary: 34C15. Secondary: 34C25,
34D10. 00footnotetext: Melnikov method, periodic orbits, completely integrable
systems..  00footnotetext:  00footnotetext: Email addresses:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
F. Crespo1, M. Uribe2, E. Martínez3
1,3GISDA, Dept. de Matemática, Facultad de Ciencias, Universidad del Bío-Bío,
Collao 1202, Casilla 5-C. Concepción, Chile
2Dept. de Matemática y Física Aplicadas, Universidad Católica de la Santísima
Concepción,
Alonso de Ribera 2850, Concepción, Chile

###### Abstract

We consider a completely integrable system of differential equations in
arbitrary dimensions whose phase space contains an open set foliated by
periodic orbits. This research analyzes the persistence and stability of the
periodic orbits under a nonlinear periodic perturbation. For this purpose, we
use the Melnikov method and Floquet theory to establish conditions for the
existence and stability of periodic orbits. Our approach considers periods of
the unperturbed orbits depending on the integrals and constant periods. In the
applications, we deal with both cases. Precisely, we study the existence of
periodic orbits in a perturbed generalized Euler system. In the degenerate
case, we analyze the existence and stability of periodic orbits for a
perturbed harmonic oscillator.

##  1 Introduction

Our research addresses the existence and stability of periodic orbits in
perturbed, completely integrable systems. Although integrable systems are rara
avis, the relevance and countless applications of systems given as perturbed
Euler and Kepler systems and the ubiquity of the harmonic oscillator justify,
in our opinion, a dedicated study for this kind of systems. Moreover, in many
applications, we are led to perturbed harmonic oscillators or pendulum systems
after normalization.

An ordinary differential system of equations is often easier to solve if we
can find one or more first integrals, usually representing conservation laws
associated with the studied phenomena. The Kepler and harmonic oscillator
systems are Liouville (Abelian) integrable [7] and maximally superintegrable
[4]. Hence, they are integrable in the broad sense introduced by
Bogoyavlenskij [2], where the author discussed Abelian and non-Abelian
integrability. In this work, we are not restricted to the Hamiltonian
formalism and consider integrability in the conventional sense. That is to
say, by completely integrable system, we mean an ordinary $N$-dimensional
system of differential equations with $(N-1)$ functionally independent first
integrals.

This work will consider the perturbations of completely integrable systems
with an open set full of periodic orbits. Then, we use the Melnikov method to
study the persistence of the mentioned periodic orbits under non-autonomous
perturbations. Precisely, we consider the following kinds of systems

$\displaystyle\begin{aligned} \dot{x}=f_{0}(x)+\epsilon
f_{1}(x,t;\epsilon),\quad x(t_{0})=x_{0},\end{aligned}$ (1)

where $x\in\mathcal{D}\subseteq\mathbb{R}^{N}$, $\mathcal{D}$ is an open set,
$f_{0}\in\mathcal{C}^{2}(\mathcal{D};\mathbb{R}^{N})$ and
$f_{1}\in\mathcal{C}^{2}(\mathcal{D}\times\mathbb{R}\times(-\delta,\delta);\mathbb{R}^{N})$
for a particular $\delta>0$. For the unperturbed case, $\epsilon=0$, the phase
space of system (1) is assumed to be full of periodic orbits. We will restrict
to the case of a periodic perturbation. That is to say, the function $f_{1}$
is $T$-periodic in the argument $t$.

Our main results provide conditions for the existence and stability of
periodic orbits in system (1). Precisely, Theorem 1 establishes conditions for
the existence of periodic solutions of system (1) by using the Implicit
Function Theorem and the Poincaré map. This theorem considers the case in
which the unperturbed orbits’ period varies with the associated integrals’
levels in system (1) and the degenerated case of a constant period. Moreover,
we provide conditions for analyzing the stability of the periodic orbits in
Theorem 2 and Corollary 1. In the applications, we establish the existence of
periodic orbits for the non-degenerate and degenerate cases, whose
perturbations contain Jacobi elliptic and trigonometric functions,
respectively. The stability is specified for the degenerate case.

The primary obstruction to the applicability of the Melnikov method is failure
to comply with the non-degeneracy conditions, leading to the consideration of
higher-order Melnikov theory. Degeneracies arise in several scenarios; one is
the case of degenerate resonances, which was previously studied in [13, 14],
where the author considered one and two degrees of freedom respectively, and
degeneracy is considered for frequencies whose first variation regarding its
corresponding action variable is null, while the second derivative is not, see
theorem 3.5 in [13], or section 2.2 in [14]. However, our treatment of
degenerate frequencies is different since we do not impose any condition on
the higher-order derivatives of the frequencies.

Oscillatory phenomena are of great importance in many fields of science and
technology. Its applications cover various issues, from engineering to
celestial mechanics or life sciences. However, our research studies an
extension of the Melnikov method as an interesting mathematical object in its
own right. This technique has proved its usefulness in very different
situations. Melnikov originally introduced it in [8] as a method to detect the
existence of transverse homoclinic points of the Poincaré map, thus allowing
the determination of chaotic dynamics through the Smale horseshoes. One year
later, this technique was also employed by Arnold [1] to prove the existence
of the phenomenon named after him as Arnold’s diffusion.  In addition, this
method is also helpful in other contexts: it allows finding sub-harmonic
periodic orbits [5, 10], which will be the objective of this paper. Moreover,
for the case of non-degenerate resonances, in [9], the author employs the
subharmonic Melnikov functions to provide sufficient conditions for real-
analytic non-integrability near periodic orbits of the unperturbed systems.

The sub-harmonic Melnikov method is usually presented in the literature for
planar systems and often deals with perturbed Hamiltonian systems, for
instance, [12, 11]. Nevertheless, due to the current development of science
and technology, real applications usually call for a high-dimensional
nonlinear system, which is frequently represented as an integrable system with
a nonlinear perturbation.

We organize our paper as follows. In Section 2, we study the structure of the
unperturbed system. Section 3 computes a standard form of the original
perturbed system based on the general solution of the unperturbed. An
approximation of the Poincaré map is given using this standard form.
Therefore, we can provide conditions for the existence and stability of
periodic orbits in Section 4. Finally, we show two applications of our results
in Section 5.

##  2 On the Geometry of the Unperturbed System

Along this section we describe the solution functions and geometric structure
of the unperturbed system ($\epsilon=0$) associated with (1).

###  2.1 Assumptions and Structure of the Unperturbed System

We start with a detailed description of the so-called unperturbed system,
which is obtained by setting $\epsilon=0$ in system (1)

$\displaystyle\begin{aligned} \dot{x}=f_{0}(x),\quad
x(t_{0})=x_{0}.\end{aligned}$ (2)

This system is an autonomous and completely integrable system of differential
equations defined in the open domain $\mathcal{D}\subseteq\mathbb{R}^{N}$.
That is to say, it is endowed with $(N-1)$ functionally independent first
integrals $\\{I_{1}(x),\ldots,I_{N-1}(x)\\}$. All these integrals are supposed
to be not locally constant and defined in a full Lebesgue measure subset
$\mathcal{D_{L}}\subset\mathcal{D}$. We impose that
$\mathcal{D}_{L}=\mathcal{D}$ in what follows. Otherwise, we may restrict the
domain of our study to the subset $\mathcal{D}_{L}$. Therefore, the following
function is a smooth $(N-1)$-submersion

$I:\mathcal{D}\rightarrow\Lambda\subset\mathbb{R}^{N-1},\quad x\rightarrow
I(x)=\left(I_{1}(x),\ldots,I_{N-1}(x)\right),$ (3)

with $\Lambda=I(\mathcal{D})$. We will use the notations $I_{x}=I(x)$ and
$\xi=(\xi_{1},...\xi_{N-1})$ for the elements in the image set $\Lambda$.
Moreover, the trajectory of every solution of system (1) is contained in
$I^{-1}(\xi)$. This description recalls the general notion of a completely
integrable system. Next, we provide the specific assumptions we shall require
from (2).

  * (i)

The pre-image $I^{-1}(\xi)\subset\mathcal{D}$ is a reunion of several disjoint
closed curves, each one is entirely traversed by a periodic solution.

  * (ii)

For each periodic orbit with initial condition $x_{0}\in\mathcal{D}$, the
period $P(x_{0})$ is a positive real number, and it is a function of type
$\mathcal{C}^{1}(\mathcal{D};\mathbb{R})$ in the argument $x_{0}$.

  * (iii)

There exists $x_{0}\in\mathcal{D}$ and non-zero $m,n\in\mathbb{N}$ such that
$m\,P(x_{0})=n\,T$, where $T$ is the period of $f_{1}$ in (1).

  * (iv)

There exists a smooth $(N-1)$-sub-manifold
$\mathcal{D}_{0}\subset\mathcal{D}$, dubbed as initial value manifold,
satisfying that each non-trivial solution of (2) intersects transversally
$\mathcal{D}_{0}$ in one and only one single point. We denote by
$x(t,x_{0}(\xi))$ the solution starting at $x_{0}(\xi)$, which is the
intersection of $I^{-1}(\xi_{1},...\xi_{N-1})$ with $\mathcal{D}_{0}$.

Prominent examples of completely integrable systems satisfying (i), (ii),
(iii) and (iv) are linear harmonic oscillators, or the Euler equations for the
free rigid body.

With the above assumptions,  the set of non-trivial solutions in the open
domain $\mathcal{D}$ is a $S^{1}$-fiber bundle space, in which the periodic
orbits of system (2) are the fibers, and the $(N-1)-$sub-manifold
$\mathcal{D}_{0}$ is the base space. Moreover, the structure of this system
imposes that the tangent space $T_{x}\mathcal{D}$ decomposes for each
$x\in\mathcal{D}$ as follows

$T_{x}\mathcal{D}=\mathcal{W}_{x}\oplus\langle u_{x}\rangle,$ (4)

where $\mathcal{W}_{x}=\langle\nabla_{x}I_{1},\ldots,\nabla_{x}I_{N-1}\rangle$
is a $(N-1)$-dimensional vector sub-space in $T_{x}\mathcal{D}$, and $u_{x}$
is a unitary vector generating $\mathcal{W}_{x}^{\bot}$.

The trajectory of a solution of system (2) with initial condition $x$ is
dubbed as the orbit through $x$, and will be denoted as $\mathcal{O}_{x}$.
Note that every orbit is contained in the following intersection

$\mathcal{O}_{x}\subseteq
I^{-1}(\xi_{1},...\xi_{N-1})=\bigcap_{j=1}^{N-1}I_{j}^{-1}(I_{x}).$ (5)

An obvious consequence of the above relation is that the vector field
$f_{0}(x)$ defining (2) belongs to the tangent bundle $T_{x}M_{j}^{x}$ of the
level hyper-manifolds
$M_{j}^{x}=\\{z\in\mathbb{R}^{N};\;I_{j}(z)=I_{j}(x)\\}$. Taking into account
that $T_{x}M_{j}^{x}=\langle\nabla_{x}I_{j}\rangle^{\bot}$, we have that

$f_{0}(x)\bot\nabla_{x}I_{j},\quad\forall j\in\\{1,..,N-1\\},$ (6)

or equivalently

$D_{x}I_{j}\cdot f_{0}(x)=0,\quad\forall j\in\\{1,..,N-1\\}.$ (7)

###  2.2 Orbits Parametrization

Let us consider a periodic orbit $\mathcal{O}_{x_{0}}$ with
$x_{0}\in\mathcal{D}_{0}$, and let $I_{0}\in\mathbb{R}^{N-1}$ be defined
through the submersion (3) as $I_{0}:=I(x_{0})$. From assumption (iv), we have
that $x_{0}$ is given by the following intersection

$x_{0}=I^{-1}(I_{0})\bigcap\mathcal{D}_{0}.$

Then, given an arbitrary $x_{0}\in\mathcal{D}_{0}$, we denote the time
parametrization of $\mathcal{O}_{x_{0}}$ as follows

$q(t,I_{0}):=x(t,x_{0}),$

where $x(t,x_{0})$ is the solution of (2) through $x_{0}$. Note that
$q(t,I_{0})$ is $P(x_{0})$-periodic. Due to the fact that $x_{0}$ and $I_{0}$
are in unique correspondence, we will express the period of
$\mathcal{O}_{x_{0}}$ as $P(I_{0}):=P(x_{0})$, which allows us to establish
the following map

$G:S^{1}\times\Lambda\subset
S^{1}\times\mathbb{R}^{N-1}\longrightarrow\mathcal{D}\subset\mathbb{R}^{N},\quad(I,\theta)\rightarrow
G(I,\theta):=q(\lambda\,\theta,I),$ (8)

where $S^{1}=\mathbb{R}/T$, $I=(I_{1},\ldots,I_{N-1})$ and $\lambda=P(I)/T$.
Note that $G$ is a diffeomorphism since it is a bijective differentiable
function. This map allows traversing all the periodic orbits in $\mathcal{D}$
at a uniformly $T$-periodic rate. Furthermore, it satisfies the following
property, which will be useful later on.

######  Lemma 1.

The partial derivatives $\\{G_{I_{1}},...,G_{I_{N-1}},G_{\theta}\\}$ generate
a basis in $\mathbb{R}^{N}$. Moreover, we have

$G_{\theta}(I,\theta)=\dfrac{P(I)}{T}f_{0}(G(I,\theta)),$ (9)

and

$DI_{i}\cdot G_{I_{j}}(I,\theta)=\delta_{j}^{i},\quad DI_{i}\cdot
G_{\theta}(I,\theta)=0,$ (10)

where $\delta_{j}^{i}$ is the Kronecker symbol and $i,j=1,...,N-1$.

###### Proof.

Since $G$ is a diffeomophism we have that the vectors
$\\{G_{\theta},G_{I_{1}},...,G_{I_{N-1}}\\}$ must be independent and therefore
a basis in $\mathbb{R}^{N}$. Now we shall prove (9) and (10). First, taking
into account that $t=\dfrac{P(I)}{T}\theta$ and differentiating $G$ with
respect to $t$ we obtain

$\dfrac{d}{dt}(G)=G_{\theta}\cdot\dot{\theta}=G_{\theta}\,\dfrac{T}{P(I)}=f_{0}(G).$

Therefore, we have $G_{\theta}(I,\theta)={P(I)}/{T}f_{0}(G(I,\theta)),$ which
proves (9). For the remaining partial derivatives we consider the identities

$I_{j}(G)=I_{j},\quad j=1,...,N-1.$

We differentiate them with respect to $I_{i}$ for $i=1,...,N-1$, which leads
to

$DI_{j}\cdot G_{I_{i}}(I,\theta)=\delta_{j}^{i}.$

Finally, $DI_{i}\cdot G_{\theta}(I,\theta)=0$ is obtained from (7). ∎

##  3 Standard Form and Melnikov Function

As it is customary, we will determine periodic orbits of the differential
system (1) by analyzing the presence of fixed points on the Poincaré map. For
this purpose, we will need our system expressed in a particular form, which we
obtain after carrying out the following suitable change of variables

$x\rightarrow(I,\theta),\quad x:=G(I,\theta),$ (11)

where $G$ is the diffeomorphic transformation given in (8). Now, we search for
the expression of system (1) in the new variables. In this regard, we consider
the solution $x=x(t,x_{0};\epsilon)$ of system (1) through
$x_{0}\in\mathcal{D}_{0}$, and differentiating with respect to $t$ we obtain

$\dot{x}=f_{0}(G)+\epsilon
f_{1}(G,t;\epsilon)=G_{I_{1}}\dot{I}_{1}+\cdots+G_{I_{N-1}}\dot{I}_{N-1}+G_{\theta}\dot{\theta}.$
(12)

Next, we derive the expressions for $\dot{\theta}$,
$\dot{I}_{1},\ldots,\dot{I}_{N-1}$, which poses the system (1) in a convenient
form for the computation of the Poincaré map. We will refer to the
differential equations $(\dot{\theta},\dot{I})$ as the standard form of system
(1).

######  Proposition 1.

After applying the change of variables given in (11), the system of
differential equation (1) reads as follows

$\displaystyle\begin{aligned}
\dot{I}_{i}=&\,\epsilon\,F_{i+1}(I,\theta,t;\epsilon),\quad i=1,...,N-1,\\\
\dot{\theta}=&\,\Omega_{I}+\epsilon\,F_{1}(I,\theta,t;\epsilon),\end{aligned}$
(13)

where $\Omega_{I}=T/P(I)$, and

$\displaystyle\begin{aligned}
F_{1}(I,\theta,t;\epsilon)=&\,\Omega_{I}\,\left[\dfrac{f_{0}(G(I,\theta))\cdot
f_{1}(G(I,\theta),t;\epsilon)}{\langle
f_{0}(G(I,\theta)),f_{0}(G(I,\theta))\rangle}-\sum_{i=1}^{N-1}\dfrac{f_{0}(G)\cdot
G_{I_{i}}(G)}{\langle f_{0},f_{0}\rangle}-DI_{i}\,f_{1}(G)\right],\\\
F_{i+1}(I,\theta,\tau;\epsilon)=&\,DI_{i}(G(I,\theta))\cdot
f_{1}(G(I,\theta),t;\epsilon),\quad i=1,...,N-1.\end{aligned}$ (14)

###### Proof.

The equation for each $\dot{I}_{i}$ is obtained by multiplying the expression
(12) by $DI_{i}$, which yields

$DI_{i}\cdot f_{0}+\epsilon DI_{i}\cdot f_{1}=DI_{i}\cdot
G_{I_{1}}\dot{I}_{1}+\cdots+DI_{i}\cdot G_{I_{N-1}}\dot{I}_{N-1}+DI_{i}\cdot
G_{\theta}\dot{\theta}.$ (15)

Taking into account Lemma 1 and the equality (7) we have that $DI_{i}\cdot
f_{0}=0$ and $DI_{i}\cdot G_{I_{j}}=\delta_{j}^{i}$. Thus, the above
expression becomes

$\epsilon DI_{i}\cdot f_{1}=\dot{I}_{i}.$ (16)

Proceeding in a similar way, we obtain the equation for $\dot{\theta}$. Now we
multiply (12) by

$\dfrac{f_{0}}{\langle f_{0},f_{0}\rangle}+\sum_{i=1}^{N-1}DI_{i}.$

Again, using Lemma 1 and (7), and after some algebraic manipulations, we get
the following relation

$1+\epsilon\dfrac{f_{0}\cdot f_{1}}{\langle
f_{0},f_{0}\rangle}=\dfrac{1}{\Omega_{I}}\dot{\theta}+\epsilon\left[\sum_{i=1}^{N-1}\dfrac{f_{0}(G)\cdot
G_{I_{i}}(G)}{\langle f_{0},f_{0}\rangle}-DI_{i}\,f_{1}(G)\right].$

Then, resolving in $\dot{\theta}$ we obtain (13). ∎

###  3.1 Fixed Points of the Poincaré Map

System (13) can be expressed as an autonomous nonlinear $(N+1)$-system by
adding an extra equation involving the independent variable as follows

$\displaystyle\begin{aligned}
\dot{I}_{i}=&\,\epsilon\,F_{i+1}(I,\theta,\tau;\epsilon),\quad i=1,...,N-1,\\\
\dot{\theta}=&\,\Omega_{I}+\epsilon\,F_{1}(I,\theta,\tau;\epsilon),\\\
\dot{\tau}=&\,1,\end{aligned}$ (17)

where the dot means $d/dt$ and $\tau\in\mathbb{R}/T$ is given by
$\tau(t)=t\;mod(T)$. In what follows, we employ the compact notation
$Z=(I,\theta,\tau)$.

######  Remark 1.

For $\epsilon=0$, the frequency $\Omega_{I}=T/P(I)$ is constant and the above
equations are easily integrated. Precisely, for the angle $\theta$ we have the
linear function $\theta(t)=\theta_{0}+\Omega_{I}t$. Thus, taking into account
that $\theta\in\mathbb{R}/T$, and assumption (iii) on section 2, we obtain
that $\theta(t)$ is periodic with minimum period equal to $nT$

$\theta(t+nT)=\theta_{0}+\Omega_{I}t+\Omega_{I}\,nT=\theta_{0}+\Omega_{I}t+m\,T=\theta_{0}+\Omega_{I}t=\theta(t).$

Since system (17) is determined by means of a $T$-periodic field in the
variables $\theta$ and $\tau$, the phase space is given by
$S^{1}\times\Lambda\times S^{1}$, with $S^{1}=\mathbb{R}/T$. Next, we are
going to construct a Poincaré map following [6]. For this aim, we define a
global cross section

$\Sigma=\\{(I,\theta,\tau)\in S^{1}\times\Lambda\times S^{1};\;\tau=0\\}.$
(18)

In view of the last equation in (17), all solutions of this autonomous
differential system cross $\Sigma$ transversely and the Poincaré map is
defined globally as

$P_{\epsilon}:\Sigma\longrightarrow\Sigma,\quad
Z_{0}=(I_{0},\theta_{0},0)\rightarrow
P_{\epsilon}(Z_{0}):=\Phi(T,Z_{0};\epsilon),$ (19)

where $\Phi(t,Z_{0};\epsilon)=(I(t),\theta(t),\tau(t))$ is the solution of
system (17) through $Z_{0}$. Note that $P_{\epsilon}$ is usually defined in a
subset of $\Sigma$, but our construction allows for a global definition in the
whole cross section $\Sigma$. Moreover, the travel time $T$ of every point in
the global cross section $\Sigma$ is the same.

Since the qualitative behaviour of system (13) can be studied by means of the
Poincaré map,  and taking into account remark 1, we focus on the fixed points
associated with the $n$-composition of $P_{\epsilon}\circ\dots\circ
P_{\epsilon}$, which will be denoted as $P_{\epsilon}^{n}$. The fixed points
of $P_{\epsilon}^{n}$ correspond to the periodic orbits of system (13).
Therefore, we will investigate the zeroes of the difference between a point in
$\Sigma$ and its image

$\widehat{\Delta}(Z_{0})=\Phi(nT,Z_{0};\epsilon)-Z_{0}.$

Because $\tau\in S^{1}=\mathbb{R}/T$ is given by $\tau(t)=t\;mod(T)$, the last
component of difference function is identically zero, hence we may simplify
$\widehat{\Delta}(Z_{0})$ by excluding the last component in the above
expression. That is to say, we consider now

$z_{0}:=(I_{0},\theta_{0}),\quad\phi(t,z_{0};\epsilon):=(I(t,z_{0};\epsilon),\theta(t,z_{0};\epsilon)),$

where $G(z_{0})\in\mathcal{D}_{0}$, and $\phi(t,z_{0};\epsilon)$ is the flow
associated with the equation (13). Thus, the existence of periodic orbits is
equivalent to analyze the zeroes of the displacement function

$\Delta(z_{0},\epsilon)=\phi(nT,z_{0};\epsilon)-z_{0}.$ (20)

The following section provides a first order approximation of $\Delta$.

###  3.2 Melnikov Function

To study the function $\Delta(z_{0},\epsilon)$ we need to compute the solution
of system (13) at $t=T$. The strategy that we will follow is based on the
computation of a suitable approximation of $\phi(t,z_{0};\epsilon)$, which is
the solution of system (13). It can be made with the aid of the parameters
dependence theorem from the fundamental theory of ordinary differential
equations. Since the vector field defining system (13) depends continuously
differentiable on the parameter $\epsilon$, the same is true for the solution
$\phi(t,z_{0};\epsilon)$. In particular, we have the following Taylor
expansion in a neighborhood of $\epsilon=0$

$\phi(t,z_{0};\epsilon)=\phi_{0}(t,z_{0})+\phi_{1}(t,z_{0})\epsilon+O(\epsilon^{2}),$

where $\phi_{0}(t,z_{0})=\phi(t,z_{0};0)$ is the solution of system (13) for
the unperturbed case $\epsilon=0$. Therefore,
$\phi_{0}(t,z_{0})=(I_{1,0}(t,z_{0}),...I_{N-1,0}(t,z_{0}),\theta_{0}(t,z_{0}))$
is given by

$\displaystyle\begin{aligned}
I_{i,0}(t,z_{0})=\,I_{i}(z_{0}),\qquad\theta_{0}(t,z_{0})=\,\theta_{0}+\Omega(I(z_{0}))\,t,\end{aligned}$
(21)

where $I_{0}(z_{0})=(I_{1}(z_{0}),...,I_{N-1}(z_{0}))\in\Lambda$ is a constant
vector. Moreover, the function

$\phi_{1}(t,z_{0})=\dfrac{\partial}{\partial\epsilon}\phi(t,z_{0};\epsilon)|_{\epsilon=0}$

can be obtained as the solution to the corresponding first variational
equation. Following [11], we write this equation as the initial value problem

$\dot{\phi}_{1}(t,z_{0})=\phi_{1w}(t,\phi_{0}(t))\phi_{1}+\phi_{1\epsilon}(t,\phi_{0}(t)),\quad\phi_{1}(0,z_{0})=0,$
(22)

being $w=(I,\theta)$, and

$\phi_{1w}(t,w)=\dfrac{\partial}{\partial
w}F(t,w;0)|_{w=\phi_{0}(t)},\quad\phi_{1\epsilon}(t,w)=\dfrac{\partial}{\partial\epsilon}F(t,w;\epsilon)|_{\epsilon=0},$

where $F(t,w;\epsilon)$ is the vector field defining (13). More specifically,
the function
$\phi_{1}(t,z_{0})=(I_{1,1}(t,z_{0}),...I_{N-1,1}(t,z_{0}),\theta_{1}(t,z_{0}))$
is the solution of system (22), which has the following expression

$\displaystyle\begin{aligned}
\dot{I}_{i,1}(t,z_{0})=&\,F_{i+1}(I_{0},\theta_{0},t;0),\\\
\dot{\theta}_{1}(t,z_{0})=&\,{\nabla_{(I,\theta)}\Omega(I_{0})}\cdot\phi_{1}(t,z_{0})+F_{1}(I_{0},\theta_{0},t;0),\end{aligned}$
(23)

or equivalently

$\displaystyle\begin{aligned}
\dot{I}_{i,1}(t,z_{0})=&\,F_{i+1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}),t;0),\\\
\dot{\theta}_{1}(t,z_{0})=&\,\sum_{i=1}^{N-1}\dfrac{\partial\Omega(I_{0})}{\partial
I_{i}}\cdot I_{i,1}(t,z_{0})+F_{1}(I_{0},\theta_{0},t;0),\\\ \end{aligned}$
(24)

where $I_{0}(t,z_{0})=(I_{1,0}(t,z_{0}),...,I_{N-1,0}(t,z_{0}))$.

Having into account that
$\phi_{1}(0,z_{0})=(I_{i,1}(0,z_{0}),\theta_{1}(0,z_{0}))=0$, we integrate the
above system to obtain a computable expresion of $\phi_{1}(t,z_{0})$

$\displaystyle\begin{aligned}
I_{i,1}(nT,z_{0})&=\,\int_{0}^{nT}DI_{i}\big{(}I_{0}(t,z_{0}),\theta_{0}(t,z_{0})\big{)}\cdot
f_{1}\big{(}I_{0}(t,z_{0}),\theta_{0}(t,z_{0}),t\big{)}\,dt,\\\
\theta_{1}(nT,z_{0})&=\,\sum_{i=1}^{N-1}\dfrac{\partial\Omega(I_{0})}{\partial
I_{i}}\int_{0}^{nT}\int_{0}^{t}DI_{i}\big{(}I_{0}(s,z_{0}),\theta_{0}(s,z_{0})\big{)}\cdot
f_{1}\big{(}I_{0}(s,z_{0}),\theta_{0}(s,z_{0}),s\big{)}\,ds\,dt\\\
&+\int_{0}^{nT}\Omega(I_{0}(t,z_{0}))\,\dfrac{f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))\cdot
f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}),t)}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}dt\\\
&-\int_{0}^{nT}\Omega(I_{0}(t,z_{0}))\,\sum_{i=1}^{N-1}\dfrac{f_{0}\cdot
G_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}\cdot
DI_{i}\,f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))dt.\end{aligned}$ (25)

We define the Melnikov function as the vector $\mathbf{M}=(M_{1},...,M_{N})$,
whose components are obtained from the above integrals

$M_{1}=I_{1,1}(nT,z_{0}),\,...,\,M_{N-1}=I_{N-1,1}(nT,z_{0}),\,M_{N}=\theta_{1}(nT,z_{0}).$

Then the first order approximation of the displacement function may be
computed from the Melnikov vector $\mathbf{M}$ as follows

$\Delta(z_{0},\epsilon)=(\epsilon M_{1},\epsilon M_{2},\ldots,\epsilon
M_{N-1},\Omega(I_{0})nT+\epsilon M_{N})+O(\epsilon^{2}).$ (26)

In the following section we prove our main result, which establishes
sufficient conditions to ensure the existence of periodic orbits. It is based
on the analysis of the Melnikov vector given above.

######  Remark 2.

Note that for the degenerate case given by ${\partial\Omega(I_{0})}/{\partial
I_{i}}=0$ for $i=1,\ldots,N-1$, the Melnikov component
$M_{N}=\theta_{1}(nT,z_{0})$ is considerably simplified. Moreover, for our
purposes, we will employ the scaled Melnikov component
$M_{N}=\theta_{1}(nT,z_{0})/\Omega(I_{0}(t,z_{0}))$, which abusing of notation
will be denoted with the same symbol. Precisely, we have

$\displaystyle\begin{aligned}
M_{N}&=\int_{0}^{nT}\dfrac{f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))\cdot
f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}),t)}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}dt\\\
&-\int_{0}^{nT}\sum_{i=1}^{N-1}\dfrac{f_{0}\cdot
G_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}-DI_{i}\,f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))dt.\end{aligned}$
(27)

##  4 Existence  and Stability of Periodic Orbits

This section present the main contribution of the present work. Namely,
conditions determining the existence and stability of periodic orbits.

###  4.1 Existence Conditions

The existence of periodic orbits is obtained through the zeroes of the
displacement function $\Delta(z_{0},\epsilon)$. However, the fact that
$\Delta(z_{0},0)$ vanishes identically implies that so does the derivative
with respect to $z_{0}$ and hence there is no use of the implicit function
theorem. To overcome this difficulty we consider the reduced displacement
function in the following theorem,  which provides the initial conditions for
periodic orbits $(I_{0}^{*},\theta_{0}^{*})\in S^{1}\times\mathbb{R}^{N-1}$ in
angle-action variables given by (11).

######  Theorem 1.

Let us assume that there exists an initial condition
$z_{0}^{*}=(I_{0}^{*},\theta_{0}^{*})$ of the unperturbed system such that
the associated orbit is $P(I_{0}^{*})$-periodic and the period satisfies the
resonance condition $m\,P(I_{0}^{*})=n\,T,$ for certain $m,n\in\mathbb{N}$.
Additionally, let us assume that one of the following condition holds

  * (a)

Non-degenerate case:

$\displaystyle\begin{aligned} &\dfrac{d\Omega_{I}}{dI}(I_{0}^{*})\neq 0,\\\
&M_{i}(I_{0}^{*},\theta_{0}^{*})=0,\quad i=1,...,N-1,\\\
&\left[\dfrac{(M_{1},...,M_{N-1},\Omega_{I})}{\partial(I_{0},\theta_{0})}\right]_{(I_{0}^{*},\theta_{0}^{*})}\neq
0.\end{aligned}$ (28)

  * (b)

Degenerate case:

$\displaystyle\begin{aligned} &\dfrac{d\Omega_{I}}{dI}(I_{0}^{*})=0,\\\
&M_{i}(I_{0}^{*},\theta_{0}^{*})=0,\quad i=1,...,N,\\\
&\left[\dfrac{(M_{1},...,M_{N})}{\partial(I_{0},\theta_{0})}\right]_{(I_{0}^{*},\theta_{0}^{*})}\neq
0.\end{aligned}$ (29)

Then, for each $0<\epsilon\ll 1$, the system (1) has a periodic solution
$x(t,x_{0}^{*}(\epsilon))$ of period near to $nT$, which initial condition
$x_{0}^{*}(\epsilon)$ is located in a neighborhood of
$x_{0}^{*}=G(z_{0}^{*}).$

###### Proof.

We have already seen that the periodic orbits of system (1) correspond to the
zeroes of the displacement function $\Delta(z_{0},\epsilon)$ given in (20).
Indeed, we can analyze the zeroes of the displacement function by means of the
following first order approximation of $\Delta(z_{0},\epsilon)$

$\Delta(z_{0},\epsilon)=(\epsilon M_{1},\epsilon M_{2},\ldots,\epsilon
M_{N-1},\Omega(I_{0})nT+\epsilon M_{N})+O(\epsilon^{2}).$

This function needs to be modified to obtain our conclusions when analyzing
the zeroes of cases (a) and (b). We will consider a modified version of the
first-order part of $\Delta$, which we will denote by $\tilde{\Delta}$ in both
cases.

  * (a)

Non-degenerate case: Zeroes of $\Delta(z_{0},\epsilon)$ are equivalent to
zeroes of the following modified displacement function

$\tilde{\Delta}(z_{0},\epsilon):=(M_{1},M_{2},\ldots,M_{N-1},\Omega(I_{0})nT)+O(\epsilon).$

The condition $M_{i}(I_{0}^{*},\theta_{0}^{*})=0$ for $i=1,...,N-1$ implies
that $\tilde{\Delta}(z_{0}^{*},0)=0$. Moreover, to ensure that
$\tilde{\Delta}(z_{0}^{*},\epsilon)=0$ in a small neighborhood of
$\epsilon=0$, we apply the implicit function theorem, which non-degeneracy
condition is given by the determinant of the Jacobian matrix denoted as
$Jac(\tilde{\Delta}(z_{0},0))$

$Det(Jac(\tilde{\Delta}(z_{0}^{*},0)))=\left[\dfrac{(M_{1},...,M_{N-1},\Omega_{I})}{\partial(I_{0},\theta_{0})}\right]_{(I_{0}^{*},\theta_{0}^{*})}\neq
0.$

Thus, we have that there exists $z_{0}^{*}(\epsilon)$ in a neighborhood of
$\epsilon=0$ satisfying $z_{0}^{*}(0)=z_{0}^{*}$ and
$\tilde{\Delta}(z_{0}^{*}(\epsilon),\epsilon)=\Delta(z_{0}^{*}(\epsilon),\epsilon)=0$,
which yields a fixed point of the $n$-Poincaré map and therefore a periodic
orbit. Finally, the initial condition of the periodic orbit is obtained
through the diffeomorphism (8) as
$x_{0}^{*}(\epsilon)=G(z_{0}^{*}(\epsilon))$.

  * (b)

Degenerate case: This case is analogous to case (a). However, the condition
${d\Omega_{I}}/{dI}(I_{0}^{*})=0$ implies that the modified displacement
function $\tilde{\Delta}$ is different from the corresponding one in case (a)

$\tilde{\Delta}(z_{0},\epsilon):=(M_{1},M_{2},\ldots,M_{N})+O(\epsilon).$

Moreover the non-degeneracy condition is given by

$Det(Jac(\tilde{\Delta}(z_{0}^{*},0)))=\left[\dfrac{(M_{1},...,M_{N})}{\partial(I_{0},\theta_{0})}\right]_{(I_{0}^{*},\theta_{0}^{*})}\neq
0.$

Again, the implicit function theorem proves the existence of the initial
condition $x_{0}^{*}(\epsilon)=G(z_{0}^{*}(\epsilon))$ leading to the periodic
orbit.

∎

######  Remark 3.

The Melnikov function $M_{N}$ arises in the degenerate case. Note that other
authors do not employ $M_{N}$ when dealing with degeneracies [13, 14]. This is
achieved by imposing several conditions on the higher-order derivatives of the
frequency $\Omega_{I}$, see theorem 3.5 in [13]. However, our treatment of
degenerate frequencies is different since we do not impose any restriction on
the frequencies.

###  4.2 Stability of Periodic Orbits

In order to study the stability of periodic orbits associated with system (1),
it is convenient to introduce the following notation.

$A(t;\epsilon)=\dfrac{\partial f_{0}}{\partial
x}(x(t,x_{0}^{*}))+\epsilon\left(\dfrac{\partial f_{10}}{\partial
x}(x(t,x_{0}^{*}),t)+\dfrac{\partial^{2}f_{0}}{\partial
x^{2}}(x(t,x_{0}^{*}))x_{1}(t,x_{0}^{*})\right),$ (30)

and $f_{10}(x(t,x_{0}^{*}),t)=f_{1}(x(t,x_{0}^{*}),t;0)$. Notice that, since
$x(t,x_{0}^{*})$ is periodic, $A(t,\epsilon)$ is a matrix with periodic
coefficients.

Next, we employ several well-known facts coming from Floquet theory to study
stability. They are collected in the following results.

######  Theorem 2.

For each small $\epsilon$, $0<\epsilon^{2}\ll 1$, we have

(i)

If the real parts of the characteristic exponents associated with
$A(t;\epsilon)$ are all negative, the periodic solution
$x(t,x_{0}^{*}(\epsilon))$ is stable. Moreover, the attraction is exponential
in a neighborhood of $x(t,x_{0}^{*}(\epsilon))$.

(ii)

If any of the real parts of the characteristic exponents associated with
$A(t;\epsilon)$ are positive, the periodic solution $x(t,x_{0}^{*}(\epsilon))$
is unstable.

###### Proof.

We consider the following change of variables, which employs the periodic
solution $x(t,x_{0}^{*}(\epsilon))$ given in theorem 1

$x=x(t,x_{0}^{*}(\epsilon))+y.$

This change of variables transform system (1) in the following maner

$\displaystyle\begin{aligned}
\dot{x}=\dot{x}(t,x_{0}^{*}(\epsilon))+\dot{y}&=f(x(t,x(t,x_{0}^{*}(\epsilon))+y,t;\epsilon)\\\
&=f(x(t,x_{0}^{*}(\epsilon)),t;\epsilon)+\dfrac{\partial f}{\partial
x}(x(t,x_{0}^{*}(\epsilon)),t;\epsilon)y+\mathcal{O}(y^{2}),\end{aligned}$

where $f(x,t;\epsilon)=f_{0}(x)+\epsilon f_{1}(x,t;\epsilon)$. Taking into
account that $x(t,x_{0}^{*}(\epsilon))$ is an actual solution of system (1),
we may simplify the above equations as follows

$\dot{y}=\dfrac{\partial f}{\partial
x}(x(t,x_{0}^{*}(\epsilon)),t;\epsilon)y+\mathcal{O}(y^{2}).$ (31)

Moreover, it happens to be that $y(t)=0$ is an equilibrium solution for system
(31), whose stability is equivalent to the stability of the periodic solution
$x(t,x_{0}^{*}(\epsilon))$. That is to say, we have proceeded by linearization
of the original system in a neighborhood of the periodic solution. Thus, the
matrix in system (31) has periodic coefficients and the theory of Floquet
applies. Therefore, we need to study the characteristic exponent associated
with the $nT$-periodic linear part in system (31), For this task, we consider
the expansion in terms of the small parameter $\epsilon$ of the periodic
solution

$x(t,x_{0}^{*}(\epsilon))=x_{0}(t,x_{0}^{*}(0))+\epsilon\,x_{1}(t,x_{0}^{*}(0))+\mathcal{O}(\epsilon^{2}),$
(32)

where $x_{0}(t,x_{0}^{*}(0))=x(t,x_{0}^{*})$ is the periodic orbit of the
unperturbed system, $x_{0}^{*}(0)=x_{0}^{*}$, and $x_{1}(t,x_{0}^{*}(0))$ is
the solution of the variational equation associated with system (31), which
now may be expressed as

$\dot{y}=\dfrac{\partial f}{\partial
x}(x(t,x_{0}^{*}),t;\epsilon)\,y+\epsilon\left[\dfrac{\partial^{2}f}{\partial
x^{2}}(x(t,x_{0}^{*}),t;\epsilon)x_{1}(t,x_{0}^{*})\right]y+\mathcal{O}(\epsilon^{2}y^{2}).$

Thus, we have $\dot{y}=A(t;\epsilon)\,y+\mathcal{O}(\epsilon^{2}y^{2}),$ and
conditions $(i)$ and $(ii)$ follow directly from Floquet theory. More
precisely, for each fixed $\epsilon$ there exists a periodic matrix
$Q(t;\epsilon)$ and a constant coefficient matrix $B_{\epsilon}$ such that
$Q(t;\epsilon)e^{tB_{\epsilon}}$ is a fundamental matrix. Moreover,
$Q(t;\epsilon)$ provides a change of variables $y=Q(t;\epsilon)z$ normalizing
the linear part of system (4.2), which becomes

$\dot{z}=B_{\epsilon}\,z+\mathcal{O}(\epsilon^{2}z^{2}),$

and the stability of the origin depends on the real parts of the eigenvectors
of $B_{\epsilon}$, which are the so called characteristic multipliers. ∎

######  Corollary 1 (Instability conditions).

Let us consider the following expression

$\lambda_{\epsilon}=\dfrac{1}{nT}\int_{t_{0}}^{t_{0}+nT}Tr\left[A(t;\epsilon)\right]dt,$
(33)

where $Tr\left[A(t;\epsilon)\right]$ is the trace of the matrix
$A(t;\epsilon)$. If $\lambda_{\epsilon}>0$, the periodic solution
$x(t,x_{0}^{*}(\epsilon))$ is unstable.

###### Proof.

Let $\lambda_{1},\ldots,\lambda_{N}$ be the characteristic multipliers
associated with $B_{\epsilon}$. The result follows from theorem 2 $(ii)$ and
the fact that

$\sum_{i=1}^{N}\lambda_{i}=\dfrac{1}{nT}\int_{t_{0}}^{t_{0}+nT}Tr\left[A(t;\epsilon)\right]dt.$

Thus, if $\lambda_{\epsilon}>0$ we have that at least one of the
characteristic multipliers has positive real part. ∎

##  5 Applications

Here we show two applications corresponding to the non-degenerate and
degenerate cases, respectively. We demonstrate the existence of periodic
orbits in both cases. Nevertheless, stability can not be guaranteed in the
non-degenerate case because of the complexity involved in the variational
equation leading to $x_{1}$ in (32).

###  5.1 Perturbed Generalized Euler System. The Non-Degenerate Case

As an application of the non-degenerate case of theorem 2, we study the
existence of periodic orbits in a perturbed generalized Euler system [3]. That
is to say, we consider the following system of differential equations

$\displaystyle\begin{aligned}
\dot{x}_{1}=&\,a_{1}\,x_{2}x_{3}+\epsilon\,p_{1}(x_{2},x_{2},x_{3},t),\\\
\dot{x}_{2}=&\,a_{2}\,x_{1}x_{3}+\epsilon\,p_{2}(x_{2},x_{2},x_{3},t),\\\
\dot{x}_{3}=&\,a_{3}\,x_{1}x_{2}+\epsilon\,p_{3}(x_{2},x_{2},x_{3},t),\\\
\end{aligned}$ (34)

where $a_{1},a_{2},a_{3}\in\mathbb{R}$ satisfy that $a_{1}a_{2}a_{3}<0$. The
perturbation introduces the dependence on the independent variable and is
given by the following periodic vector

$\displaystyle\begin{aligned}
p_{1}(x_{1},x_{2},x_{3},t)=&\,x_{1}+\alpha_{1}\,\text{cn}(\omega_{1}t,k)+\beta_{1}\,\text{sn}(\omega_{2}t,k),\\\
p_{2}(x_{1},x_{2},x_{3},t)=&\,x_{2}+\alpha_{2}\,\text{cn}(\omega_{1}t,k)+\beta_{2}\,\text{sn}(\omega_{2}t,k),\\\
p_{3}(x_{1},x_{2},x_{3},t)=&\,{x_{3}}+{a}\,{\text{nd}(\omega_{3}t,k)}+b\,\dfrac{x_{1}}{x_{3}}\,\text{cn}(\omega_{2}t,k),+c\,\dfrac{x_{2}}{x_{3}}\,\text{sn}(\omega_{2}t,k),\\\
\end{aligned}$ (35)

where $\alpha_{1},\beta_{1},\alpha_{2},\beta_{2},a,b,c\in\mathbb{R}$.
Moreover, $\text{cn}(t,k)$, $\text{sn}(t,k)$ and $\text{dn}(t,k)$ denote the
Jacobi elliptic functions with $k$ the elliptic modulus and using Glaisher
notation $\text{nd}(t,k)$ is the inverse of $\text{dn}(t,k)$. Note that system
(34) depends on ten parameters and involves elliptic functions.

The above example does not resemble any specific application. Instead, we
chose a non-trivial perturbation involving elliptic functions to illustrate
the previous theory. By this choice, we show that our development tackles a
wide range of perturbations, and even for the case of elliptic perturbations,
technical difficulties associated with the integrals in the Melnikov functions
are affordable. A similar and very interesting example from the applied point
of view was studied in [15], a quad-rotor helicopter model involving a
trigonometric perturbation. In a forthcoming research, we will focus on
particular physically meaningful applications.

####  5.1.1 A Précis on the Generalized Euler System

In order to have a good underestanding of the dynamics of the unperturbed
system, this section gives a brief description of the main facts about the
case $\epsilon=0$ in system (34)

$\displaystyle\begin{aligned} \dot{x}_{1}=&a_{1}\,x_{2}x_{3},\\\
\dot{x}_{2}=&a_{2}\,x_{1}x_{3},\\\ \dot{x}_{3}=&a_{3}\,x_{1}x_{2}.\\\
\end{aligned}$ (36)

This system is dubbed as Generalized Euler system in [3], and it is studied
for arbitrary parameters $a_{1},a_{2},a_{3}\in\mathbb{R}$. However,
interchanging variables and equations, if needed, and by choosing a suitable
orientation of the independent variable, if it is needed, system (36) may be
arranged in such a way that the parameter $a_{1}>0$, $a_{3}\geq 0$ and
$a_{2}\in\mathbb{R}$. In what follows, we will assume that our system is
expressed in this way.

It is well known that a dynamical system fitting (36) are Euler equations,
i.e., the reduced system associated with the motion of a rigid body with a
fixed point in the absence of external forces. Then, the variables $x_{i}$ are
the three components of the angular momentum, and parameters $a_{i}$ depend on
the principal moments of inertia of the rigid body.

The sign of the parameters determine two possible cases of bounded or
unbounded orbits. In this application the solutions of the unperturbed system
must be periodic. For that reason we restrict to the case $a_{1}a_{2}a_{3}<0$,
which implies that all the solutions are bounded. Next, we will impose
additional conditions to ensure that they are indeed periodic.

A simple computation shows that system (36) is endowed with the following
polynomial integrals

$h_{1}=a_{2}x_{3}^{2}-a_{3}x^{2}_{2},\quad
h_{2}=a_{3}x_{1}^{2}-a_{1}x^{2}_{3},\quad
h_{3}=a_{1}x_{2}^{2}-a_{2}x^{2}_{1}.$ (37)

For the case $a_{1}a_{2}a_{3}<0$, the integrals $h_{1}$ and $h_{3}$ are two
elliptic cylinders and a $h_{2}$ is a hyperbolic one, whose intersection gives
the trajectories of the solutions for system (36).

In [3], it is proven that given initial conditions such that $h_{1}\neq 0$,
$h_{2}\neq 0$ and $h_{3}\neq 0$, the solutions trajectories of system (36) are
not contained but equal to each of the connected component of the intersection
of the integrals. For the case $a_{2}<0$ we obtain periodic and bounded
solutions, and $a_{2}>0$ leads to the non-bounded ones. If any of the
parameters vanishes, we have the following cases; $a_{j}=0$ and $a_{i}a_{k}>0$
leads to unbounded solutions, while $a_{j}=0$, $a_{i}a_{k}<0$ to the bounded
ones. The case in which two parameters vanish corresponds to straight
trajectories.

Note that the integral $h_{2}$ has a key role separating two different types
of dynamics. Precisely for $h_{2}<0$ we have oscillations around the
$x_{3}$-axis, for $h_{2}>0$ the oscillations are around $x_{1}$-axis, and
$h_{2}=0$ separates both dynamics by herteroclinic orbits.

Finally, we provide explicitly the solution of system (36) for the bounded
case. Moreover, our formulas are restricted for the case $a_{2}h_{a}>0$.
Namely, we have that the solutions of the bounded extended Euler system with
initial conditions

$x_{1}(s_{0})=x_{1}^{0}=\sqrt{\frac{h_{3}}{a_{2}}},\quad
x_{2}(s_{0})=x_{2}^{0}=0,\quad
x_{3}(s_{0})=x_{3}^{0}=\sqrt{\frac{h_{1}}{a_{2}}},$

may be expressed as

$\displaystyle\begin{aligned}
&x_{1}(t+t_{0})=\phantom{-}\sqrt{\Big{|}\dfrac{h_{3}}{a_{2}}\Big{|}}{\text{c}n}(\mu\,t,k),\\\
&x_{2}(t+t_{0})=-\sqrt{\Big{|}\dfrac{h_{3}}{a_{1}}\Big{|}}{\text{s}n}(\mu\,t,k),\\\
&x_{3}(t+t_{0})=\phantom{-}\sqrt{\Big{|}\dfrac{h_{1}}{a_{2}}\Big{|}}{\text{d}n}(\mu\,t,k),\end{aligned}$
(38)

where the elliptic modulus is given by

$k=\sqrt{\Big{|}\dfrac{a_{3}h_{3}}{a_{1}h_{1}}\Big{|}},\qquad\mu=\sqrt{|a_{1}h_{1}|}.$

Furthermore, $x_{3}(s)$ and $x_{2}(s)$ are $4\mathcal{P}$ periodic and
$x_{1}(s)$ is $2\mathcal{P}$ periodic if $a_{2}h_{2}>0$. The period may be
calculated through $\mathcal{P}$ as the integral given by

$\mathcal{P}(a_{1},h_{1},a_{3},h_{3})=\int_{0}^{1}\dfrac{\sqrt{|a_{1}/h_{1}|}}{\sqrt{(1-z^{2})(1-k^{2}z^{2})}}\,dz.$
(39)

Notice that the above solutions (LABEL:eq:BoundedSolutionsGeneral) are
periodic orbits oscillating around the $x_{3}$-axis.

####  5.1.2 Perturbed Orbits Throught Melnikov Method

Now we will show the existence of periodic orbits associated with the
perturbed generalized Euler system given in (34). For this purpose, we will
fix the value of the parameters to

$a_{1}=a_{3}=1,\;a_{2}=-2,\quad\alpha_{1}=-\dfrac{3}{2\sqrt{2}},\;\beta_{1}=2.89972,\quad\alpha_{2}=0,\;\beta_{2}=\dfrac{3}{2},$
$a=-\dfrac{3}{\sqrt{2}},\;b=0,\;c=-\dfrac{3}{4}\quad\omega_{1}=1,\;\omega_{2}=2,\;\omega_{3}=3,\;k=\dfrac{1}{2},$

which yields

$\displaystyle\begin{aligned}
\dot{x}_{1}=&\,x_{2}x_{3}+\epsilon\,\left(x_{1}-\dfrac{3}{2\sqrt{2}}\text{cn}(t,\frac{1}{2})+2.89972\,\text{sn}(2t,\frac{1}{2})\right),\\\
\dot{x}_{2}=&-2\,x_{1}x_{3}+\epsilon\,\left(x_{2}+\dfrac{3}{2}\text{sn}(2t,\frac{1}{2})\right),\\\
\dot{x}_{3}=&\,x_{1}x_{2}+\epsilon\,\left({x_{3}}-\dfrac{3}{\sqrt{2}}\,\text{nd}(3t,\frac{1}{2})-\dfrac{3}{4}\dfrac{x_{2}}{x_{3}}\,\text{sn}(2t,\frac{1}{2})\right).\\\
\end{aligned}$ (40)

Figure 1: Periodic orbit simulation. Figure 1 and Figure 1 shows the periodic
orbit for $\epsilon=0.001$ and $\epsilon=0.002$. In Figure 1, we plot the
evolution of $Norm[(x_{1}(t),x_{2}(t),x_{3}(t))-(x_{1}(0),x_{2}(0),x_{3}(0))]$
for eleven cycles and $\epsilon=0.001$.

The unperturbed system $\epsilon=0$ satisfies the hypothesis (i), (ii) and
(iii) assumed in Section 2. Precisely, we define the open set
$\mathcal{D}\subset\mathbb{R}^{3}$ by

$\mathcal{D}=\\{(x_{1},x_{2},x_{3})\in\mathbb{R}^{3};\;a_{2}h_{2}>0,\,x_{3}>0\\},$

which is made only of periodic orbits. The so called initial value manifold
$\mathcal{D}_{0}$ is the intersection of $\mathcal{D}$ with the plane
$x_{2}=0$, and the vector of integrals $I=(I_{1},I_{2})$ is given by

$I_{1}=\sqrt{|h_{1}/a_{1}|},\quad I_{2}=\sqrt{|h_{3}/a_{3}|}.$

The period $P(x_{0})=P(I)$ of any solution starting in $\mathcal{D}$ is
computed by (39) as

${P}(I_{1},I_{2})=\dfrac{4}{\sqrt{I_{1}}}\int_{0}^{1}\dfrac{1}{\sqrt{(1-z^{2})(1-\frac{I_{2}}{I_{1}}^{2}z^{2})}}\,dz.$
(41)

Therefore, ${P}(I_{1},I_{2})$ is a smooth function. It is equal to the period
of the perturbation for any initial condition in the curve obtained by the
intersection $I_{1}=2$ and $I_{2}=1$. Now, using the solution of the
unperturbed system (LABEL:eq:BoundedSolutionsGeneral) and recalling that we
denote $T$ for the period of the perturbation, we have that the change of
variables $(x_{1},x_{2},x_{3})\rightarrow(I_{1},I_{2},\theta)$ given in (11)
is determined by

$\displaystyle\begin{aligned}
x_{1}(I_{1},I_{2},\theta)=&I_{2}\,\text{cn}(\sqrt{2}\,I_{1}\lambda\,\theta,\frac{I_{2}}{I_{1}}),\\\
x_{2}(I_{1},I_{2},\theta)=&-\sqrt{2}\,I_{2}\,\text{sn}(\sqrt{2}I_{1}\lambda\,\theta,\frac{I_{2}}{I_{1}}),\\\
x_{3}(I_{1},I_{2},\theta)=&I_{1}\,\text{dn}(\sqrt{2}\,I_{1}\lambda\,\theta,\frac{I_{2}}{I_{1}}),\end{aligned}$
(42)

where $\lambda=P(I)/T$.

Taking into account that $\Omega_{I}=1/\lambda=T/P(I_{1},I_{2})$ correspond
with the non-degenerate case, we compute the Melnikov vector $(M_{1},M_{2})$
by means of equations (25)

$\displaystyle\begin{aligned}
M_{1}=&\,\dfrac{1}{I_{1}}\int_{0}^{T}x_{2}^{2}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)+2x_{3}^{2}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)+\dfrac{-3\sqrt{2}\,x_{3}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)}{\text{dn}(3t,\frac{1}{2})}\,dt,\\\
M_{2}=&\,\dfrac{2}{I_{2}}\int_{0}^{T}{x_{1}}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)\left(\frac{-3}{2\sqrt{2}}\text{cn}(t,{1}/{2})+2.89972\,\text{sn}(2t,{1}/{2})+{x_{1}}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)\right)\\\
&\phantom{\,\dfrac{2}{I_{2}}\int_{0}^{T}}+{x_{2}}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)\left(\frac{3}{2}\,\text{sn}(2t,{1}/{2})+{x_{2}}(I_{1},I_{2},\Omega_{I}t+\theta_{0},)\right)\,dt.\end{aligned}$
(43)

After some cumbersome computations, that involves the integration with several
Jacobi elliptic functions, we obtain the following zero of the Melnikov vector

$\theta_{0}^{*}=0,\quad I_{1}^{*}=\dfrac{3}{\sqrt{2}},\quad
I_{2}^{*}=\dfrac{3}{2\sqrt{2}},$

with the following degeneration matrix

$\begin{bmatrix}\partial\Omega_{I}/\partial\theta_{0}&\partial\Omega_{I}/\partial
I_{1}&\partial\Omega_{I}/\partial I_{2}\\\ \partial
M_{1}/\partial\theta_{0}&\partial M_{1}/\partial I_{1}&\partial M_{1}/\partial
I_{2}\\\ \partial M_{2}/\partial\theta_{0}&\partial M_{2}/\partial
I_{1}&\partial M_{2}/\partial I_{2}\end{bmatrix}=-5.16.$ (44)

Thus, we can ensure that system (40) has a periodic orbit for small $\epsilon$
in the neighborhood of the initial condition

$(x_{1}(I_{1}^{*},I_{2}^{*},\theta_{0}^{*}),x_{2}(I_{1}^{*},I_{2}^{*},\theta_{0}^{*}),x_{3}(I_{1}^{*},I_{2}^{*},\theta_{0}^{*}))=(1.06066,0,2.12132).$

In Figure 1 we plot the numerical evaluation of the solution of system (40)
for the above initial condition. It illustrates the periodic orbit for two
small values of $\epsilon$ and also shows the behavior for several cycles of
the oscillation.

###  5.2 Existence and Stability in the Degenerate Case

Let us consider the following system of differential equations

$\displaystyle\begin{aligned}
\dot{x}_{1}=&\,\omega\,y_{1}+\epsilon\,a\,\cos(t),\\\
\dot{y}_{1}=&\,-\omega\,x_{1}+\epsilon
b\left(x_{1}\left(x_{2}^{2}+y_{2}^{2}\right)\cos(t)\right),\\\
\dot{x}_{2}=&\,\omega\,y_{2}+\epsilon\,c\,x_{1}\left(1-\left(x_{2}^{2}+y_{2}^{2}\right)\right),\\\
\dot{y}_{2}=&\,-\omega\,x_{1}+\epsilon\,d\,x_{1}\left(1-\left(x_{2}^{2}+y_{2}^{2}\right)\right),\\\
\end{aligned}$ (45)

where $\omega=1$ and $\\{a,b,c,d\\}=\\{601/(6\sqrt{5}),0,-1/1080,1/2280\\}$.
This system is endowed with the following quadratic homogeneous polynomials
integrals

$I_{1}=x_{1}^{2}+y_{1}^{2},\quad I_{2}=x_{2}^{2}+y_{2}^{2},\quad
I_{3}=x_{1}x_{2}+y_{1}y_{2}.$

For $\epsilon=0$, the general solution of the above system reads as follows

$\displaystyle\begin{aligned}
x_{1}&=\sqrt{I_{1}}\cos(\theta_{0}+\omega\,t),\\\
y_{1}&=\sqrt{I_{1}}\sin(\theta_{0}+\omega\,t)\\\
x_{2}&=\frac{I_{3}\cos(\theta_{0}+\omega\,t)-\sqrt{I_{1}I_{2}-I_{3}^{2}}\sin(\theta_{0}+\omega\,t)}{\sqrt{I_{1}}},\\\
y_{2}&=\frac{I_{3}\sin(\theta_{0}+\omega\,t)+\sqrt{I_{1}I_{2}-I_{3}^{2}}\cos(\theta_{0}+\omega\,t)}{\sqrt{I_{1}}}.\end{aligned}$
(46)

These formulas provide us with the change of variables $x=G(I,\theta)$ given
in (11), where $x=\\{x_{1},y_{1},x_{2},y_{2}\\}$. Moreover, in our case
$nT=2\pi$, and the Melnikov vector is given by
$\\{M_{1},M_{2},M_{3},M_{4}\\}$, where

$\displaystyle\begin{aligned}
M_{1}=\int_{0}^{nT}&2b\,x_{1}\cos(t)+2y_{1}\left(x_{1}\left(x_{2}^{2}+y_{2}^{2}\right)+a\cos(t)\right)dt,\\\
M_{2}=\int_{0}^{nT}&2c\,x_{1}x_{2}\left(-x_{2}^{2}-y_{2}^{2}+1\right)+2d\,x_{1}y_{2}\left(-x_{1}^{2}-y_{1}^{2}+1\right)dt\\\
M_{3}=\int_{0}^{nT}&y_{2}\left(a\cos(t)+x_{1}\left(x_{2}^{2}+y_{2}^{2}\right)\right)+b\,x_{2}\cos(t).\end{aligned}$
(47)

The expression for $M_{4}$ is computed by taking into account remark 2

$\displaystyle\begin{aligned}
M_{4}&=\int_{0}^{nT}\dfrac{f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))\cdot
f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}),t)}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}dt\\\
&-\int_{0}^{nT}\sum_{i=1}^{3}\dfrac{f_{0}\cdot
G_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))}{|f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))|^{2}}-DI_{i}\,f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))dt,\end{aligned}$
(48)

where $f_{0}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))$ and
$f_{1}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))$ are obtained by plugin formulas
(46) into

$\displaystyle\begin{aligned} f_{0}&=\omega\\{y_{1},-x_{1},y_{2},-x_{2}\\}\\\
f_{1}&=\omega\\{a\cos(t),0,c\,x_{1}\left(1-\left(x_{2}^{2}+y_{2}^{2}\right)\right),d\,x_{1}\left(1-\left(x_{2}^{2}+y_{2}^{2}\right)\right)\\}.\end{aligned}$
(49)

Moreover, we compute $G_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))$ and
$D_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))$ as

$\displaystyle\begin{aligned}
G_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))&=\dfrac{\partial}{\partial
I_{i}}\\{x_{1}(I,\theta),y_{1}(I,\theta),x_{2}(I,\theta),y_{2}(I,\theta)\\},\end{aligned}$
(50)

while $D_{I_{i}}(I_{0}(t,z_{0}),\theta_{0}(t,z_{0}))$ is obtained as the
differential of $I_{i}(x_{1},y_{1},x_{2},y_{2})$, and then plugin formulas
(46). We found the following zero of the Melnikov function
$(I_{0}^{*},\theta_{0}^{*})=(I_{1}^{*},I_{2}^{*},I_{3}^{*},\theta_{0}^{*})$,
where $I_{1}^{*}=5,\quad I_{2}^{*}=2,\quad
I_{3}^{*}=1,\quad\theta_{0}^{*}=\pi$, which satisfies that

$\begin{bmatrix}\partial M_{1}/\partial\theta_{0}&\partial M_{1}/\partial
I_{1}&\partial M_{1}/\partial I_{2}&\partial M_{1}/\partial I_{3}\\\ \partial
M_{2}/\partial\theta_{0}&\partial M_{2}/\partial I_{1}&\partial M_{2}/\partial
I_{2}&\partial M_{2}/\partial I_{3}\\\ \partial
M_{3}/\partial\theta_{0}&\partial M_{3}/\partial I_{1}&\partial M_{3}/\partial
I_{2}&\partial M_{3}/\partial I_{3}\\\ \partial
M_{4}/\partial\theta_{0}&\partial M_{4}/\partial I_{1}&\partial M_{4}/\partial
I_{2}&\partial M_{4}/\partial I_{3}\end{bmatrix}=49.9940.$ (51)

Thus, for each $\epsilon$, there is a periodic orbit of system (45) with
initial condition near to

$\left(x_{1}^{*}(I_{0}^{*},\theta_{0}^{*}),y_{1}^{*}(I_{0}^{*},\theta_{0}^{*}),x_{2}^{*}(I_{0}^{*},\theta_{0}^{*}),y_{2}^{*}(I_{0}^{*},\theta_{0}^{*})\right)=\left(-2\sqrt{5},-\sqrt{5},0,-\sqrt{5}\right).$

Figure 2: Evolution of
$Norm[(x_{1}(t),y_{1}(t),x_{2}(t),y_{2}(t))-(x_{1}(0),y_{1}(0),x_{2}(0),y_{2}(0))]$
for 72 periods. In Figure 2 we consider $\epsilon=0.1$ and Figure 2 is for
$\epsilon=-0.1$.

Figure 3: Periodic orbit simulation: we plot the projection
$(x_{1}(t),x_{2}(t),y_{1}(t))$ for 72 periods. In Figure 3 we consider
$\epsilon=0.1$ and Figure 3 is for $\epsilon=-0.1$.

Moreover, we assess the stability by means of corollary 1, which is in this
case yields

$\dfrac{1}{2\pi}\int_{t_{0}}^{t_{0}+2\pi}Tr\left[A(t;\epsilon)\right]dt=0.0054\,\epsilon.$

Therefore, we conclude that the periodic orbits are unstable for $\epsilon>0$.
Next, we illustrate our conclusions with a numeric simulation for $\epsilon>0$
and $\epsilon<0$. Figure 2 shows the evolution of
$Norm[(x_{1}(t),y_{1}(t),x_{2}(t),y_{2}(t))-(x_{1}(0),y_{1}(0),x_{2}(0),y_{2}(0))]$
for 72 periods. In the case $\epsilon<0$ and for an initial condition near
$x^{*}(I_{0}^{*},\theta_{0}^{*})$ we have a trajectory with a very regular
behavior. However, for $\epsilon>0$, the solution begin to tumble after 8
periods and finally scapes for more that 72. This behavior is confirmed in
Figure 3, where, for the same initial condition, we plot a close to periodic
trajectory for $\epsilon<0$ and 72 periods. The case $\epsilon>0$ shows an
orbit which eventually moves away from the periodic orbit.

## Acknowledgements

Support from the research projects GISDA 196108 GI/C UBB is acknowledged. The
author F.C. acknowledges support from 2020157 IF/R UBB. The second author was
partially suported by Proyecto DIREG 14/2022, Dirección de Investigación de la
Universidad Católica de la Santísima Concepción, Chile and Proyecto Ingeniería
2030 (ING222010004).

## References

  * [1] V. I. Arnold. Instability of dynamical systems with many degrees of freedom. Dokl. Akad. Nauk SSSR, 156(1):9–12, 1964.  doi.org/10.1007/978-3-642-01742-1_26.
  * [2] O. I. Bogoyavlenskij. Extended integrability and bi-Hamiltonian systems. Comm. Math. Phys., 196:19–51, 1998.  doi.org/10.1007/s002200050412.
  * [3] F. Crespo and S. Ferrer. On the extended Euler system and the Jacobi and Weierstrass elliptic functions. Journal of Geometric Mechanics, 7(2):151–168, 2015.  doi.org/10.3934/jgm.2015.7.151.
  * [4] F. Fassò. Superintegrable Hamiltonian systems: geometry and perturbations. Acta Appl. Math., 87:93–121, 2005.  doi.org/10.1007/s10440-005-1139-8.
  * [5] B. D. Greenspan and P. Holmes. Homoclinic orbits, subharmonics and global bifurcations in forced oscillations. Nonlinear Dynamics and Turbulence, G. Barenblatt, B. Iooss, and D. D. Joseph, eds., Pitman, London, pp. 172-214, 1983.
  * [6] J. Guckenheimer and P. Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer, New York, 1983.
  * [7] J. Liouville. Note sur l’intégration des équations différentielles de la dynamique. J. Math. Pures Appl., 20:137–138, 1855.
  * [8] V. K. Mel’nikov. On the stability of a center for time-periodic perturbations. Tr. Mosk. Mat. Obs., 12:3–5, 1963.
  * [9] S. Motonaga. Subharmonic Melnikov functions and nonintegrability for autonomous and non-autonomous perturbations of single-degree-of-freedom Hamiltonian systems near periodic orbits. Physica D: Nonlinear Phenomena, 460:134088, 2024.  doi.org/10.1016/j.physd.2024.134088.
  * [10] M. Sun, W. Zhang, and J.E. Chen. Subharmonic Melnikov method of six-dimensional nonlinear systems and application to a laminated composite piezoelectric rectangular plate. Nonlinear Dynamics, 75(1):289–310, 2014.  doi.org/10.1007/s11071-013-1066-6.
  * [11] G. Teschl. Ordinary differential equations and dynamical systems. American Mathematical Society, 2012.
  * [12] S. Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer-Verlag New York, 2003.
  * [13] K. Yagasaki, The Melnikov theory for subharmonics and their bifurcations in forced oscillations, SIAM J. Appl. Math. 56(6), 1720–1765, 1996.  doi.org/10.1137/S0036139995281317.
  * [14] K. Yagasaki, Degenerate resonances in forced oscillators, Discrete and Continuous Dynamical Systems Ser. B. 3(3), 423–438, 2003.  doi.org/10.3934/dcdsb.2003.3.423.
  * [15] K. Yagasaki, Heteroclinic Transition Motions in Periodic Perturbations of Conservative Systems with an Application to Forced Rigid Body Dynamics, Regul. Chaot. Dyn. 23, 438–457, 2018.  doi.org/10.1134/S1560354718040056.