# The slow divergence integral and torus knots

Renato Huzak  Hasselt University, Campus Diepenbeek, Agoralaan Gebouw D, 3590
Diepenbeek, Belgium<EMAIL_ADDRESS>and  Hildeberto Jardón-
Kojakhmetov  University of Groningen, Faculty of Science and Engineering,
Dynamical Systems, Geometry & Mathematical Physics — Bernoulli Institute,
Nijenborgh 9, 9747 AG, Groningen, The Netherlands<EMAIL_ADDRESS>
###### Abstract.

The goal of this paper is to study global dynamics of $C^{\infty}$-smooth
slow-fast systems on the $2$-torus of class $C^{\infty}$ using geometric
singular perturbation theory and the notion of slow divergence integral. Given
any $m\in\mathbb{N}$ and two relatively prime integers $k$ and $l$, we show
that there exists a slow-fast system $Y_{\epsilon}$ on the $2$-torus that has
a $2m$-link of type $(k,l)$, i.e. a (disjoint finite) union of $2m$ slow-fast
limit cycles each of $(k,l)$-torus knot type, for all small $\epsilon>0$. The
$(k,l)$-torus knot turns around the $2$-torus $k$ times meridionally and $l$
times longitudinally. There are exactly $m$ repelling canard limit cycles and
$m$ attracting non-canard limit cycles.

##  1\. Introduction

Singularly perturbed systems on the $2$-torus $\mathbb{T}^{2}$ have been
studied in [7] and [9, 10, 11, 12]. In [7] the authors constructed a slow-fast
system on $\mathbb{T}^{2}$ (depending only on a singular parameter $\epsilon$)
with the following property: there is a sequence of $\epsilon$-intervals
accumulating at $0$ such that the system has exactly $2$ limit cycles (one is
a stable canard and the other one is an unstable canard) for each $\epsilon>0$
from these intervals. A limit cycle is called a canard111There exists a more
classical definition of canard: a canard limit cycle passes near both
attracting and repelling portions of the critical manifold. limit cycle if it
contains a part passing near repelling portions of the critical or slow curve.
In [9, 10, 11, 12], I. V. Schurov generalized the results of [7] (the main
focus has been directed towards the existence of (attracting) canard cycles on
$\mathbb{T}^{2}$ under more general conditions). The main tool is the Poincaré
map from $\mathbb{S}^{1}$ to itself. If the rotation number of the Poincaré
map is an integer and the critical manifold is connected, then the number of
canard limit cycles is bounded by the number of fold points of the critical
manifold ([10]).

The papers mentioned above mainly deal with “unknotted” limit cycles on the
$2$-torus (the limit cycles make one pass along the slow direction, i.e. the
case of integer rotation number) and the connected critical curve with jump
contact points does not turn around the $2$-torus along the slow direction and
the fast direction. An exception is [9] where a result was proved for non-
integer rotation number (more precisely, canard limit cycles make two passes
along the slow direction). The main purpose of our paper is to show the
existence of slow-fast systems, with one parameter $\epsilon$, with an
arbitrary finite number of repelling canard cycles on $\mathbb{T}^{2}$ that
make $k$ passes along the slow direction and $l$ passes along the fast
direction ($k$ and $l$ are relatively prime). Such limit cycles are called
$(k,l)$-torus knots and occur for each small $\epsilon>0$ (for more details
about the torus knots see the rest of this section and Appendix A). They are
generated by an unconnected normally hyperbolic critical curve (each component
is a $(k,l)$-torus knot). In order to prove this and to find the global
dynamics on $\mathbb{T}^{2}$, we will use a fixed point theorem for segments,
the notion of slow divergence integral defined along the critical curve and a
generalization of Poincaré-Bendixon theorem to $\mathbb{T}^{2}$ due to
Schwartz (see [13]). In our slow-fast setting, it is more convenient to not
use the Poincaré map (see Remark 1).

To start fixing ideas, let us consider a slow-fast system

(1) $X_{\epsilon}:\left\\{\begin{array}[]{rcl}\dot{x}&=&\sin{(y-x)}\\\
\dot{y}&=&\epsilon\end{array}\right.$

where $(x,y)\in\mathbb{T}^{2}=\mathbb{R}^{2}\setminus(2\pi\mathbb{Z}^{2})$ and
$\epsilon\geq 0$ is a (small) singular parameter. The vector field
$X_{\epsilon}$ is $2\pi$-periodic in both variables and we restrict our
attention to the dynamics of $X_{\epsilon}$, with $\epsilon\geq 0$, on the
two-dimensional torus $\mathbb{T}^{2}$ (we keep $x,y$ in $[0,2\pi[$ and glue
together the opposite segments $x=0$ and $x=2\pi$, and $y=0$ and $y=2\pi$). In
the limit $\epsilon=0$, system (1) has horizontal fast orbits and two disjoint
simple closed curves of singularities given by $C_{-}=\\{y=x\\}$ and
$C_{+}=\\{y=x+\pi\\}$. These two closed curves pass through $\mathbb{T}^{2}$
horizontally and vertically only once (see Figure 1). All the singularities on
$C_{-}$ (resp. $C_{+}$) are normally attracting (resp. normally repelling).
When $\epsilon>0$, these singularities disappear and the dynamics near
$C_{\pm}$ are given by the reduced flow $x^{\prime}=1$ (or equivalently,
$y^{\prime}=1$).

$0$$C_{+}$$C_{-}$$C_{+}$$2\pi$$2\pi$

Figure 1. Dynamics of $X_{0}$.

We are interested in the dynamics of the regular system $X_{\epsilon}$ on
$\mathbb{T}^{2}$ for each small and positive parameter $\epsilon$. Any orbit
$\mathcal{O}_{\epsilon}$ of $X_{\epsilon}$ with initial point located away
from $C_{\pm}$ is attracted to $C_{-}$ and stays close to $C_{-}$ forever.
Such orbit $\mathcal{O}_{\epsilon}$ cannot therefore be closed. It is clear
now that closed orbits of $X_{\epsilon}$ may appear only in a tubular
neighborhood of $C_{-}$ or $C_{+}$. Notice that $C_{-}$ and $C_{+}$ are limit
periodic sets at level $\epsilon=0$ without fast segments ($C_{-}$ consists of
one attracting slow part while $C_{+}$ has one repelling slow part).

We show that $C_{-}$ (resp. $C_{+}$) generates precisely one limit cycle which
is hyperbolic and attracting (resp. repelling and canard), for each small and
positive $\epsilon$. This result will be true in an $\epsilon$-uniform
neighborhood of $C_{-}$ (resp. $C_{+}$). It means that the neighborhood does
not shrink to $C_{\pm}$ as $\epsilon\to 0$. The $\omega$-limit set (resp. the
$\alpha$-limit set) of all other orbits (different from the two limit cycles)
is the attracting (resp. repelling) limit cycle Hausdorff close to $C_{-}$
(resp. $C_{+}$).

This global result will be proved not only for system (1), but more general
$C^{\infty}$-smooth slow-fast systems on $\mathbb{T}^{2}$ as well. Instead of
two critical curves one of which is normally attracting and the other one is
normally repelling, we can have $m$ normally attracting closed critical curves
and $m$ normally repelling closed critical curves where $m\in\mathbb{N}$ is
arbitrary and fixed. More precisely, let us look at the following
generalization of system (1):

(2) $Y_{\epsilon}:\left\\{\begin{array}[]{rcl}\dot{x}&=&\sin{m(ly-kx)}\\\
\dot{y}&=&\epsilon\end{array}\right.$

$0$$\frac{\pi}{ml}$$\frac{\pi}{mk}$$\frac{2\pi}{mk}$$2\pi$$2\pi$$\frac{(2ml-1)\pi}{ml}$

Figure 2. Dynamics of $Y_{0}$.

where $m\in\mathbb{N}$, $k\in\mathbb{N}$, $l$ is a non-negative integer and
the pair $(k,l)$ is relatively prime. It is not difficult to see that $Y_{0}$
has $2m$ disjoint simple closed curves of singularities on $\mathbb{T}^{2}$
and that each closed curve wraps $k$ times vertically around $\mathbb{T}^{2}$
and $l$ times horizontally. This means that each closed curve crosses the
horizontal interval $[0,2\pi[\times\left\\{0\right\\}$ exactly $k$ times and
the vertical interval $\left\\{0\right\\}\times[0,2\pi[$ $l$ times (see Figure
2). Therefore, we have $m$ normally attracting curves denoted by
$C_{-}^{1},\dots,C_{-}^{m}$ and $m$ normally repelling curves denoted by
$C_{+}^{1},\dots,C_{+}^{m}$ (note that $k,m>0$). We shall show that the slow-
fast system (2) has exactly $m$ hyperbolically attracting (non-canard) limit
cycles and $m$ hyperbolically repelling (canard) limit cycles for each
$\epsilon>0$ and $\epsilon\sim 0$.

######  Theorem 1.1.

There exist $\epsilon_{0}>0$ and tubular neighborhoods $\mathcal{U}_{\pm}^{i}$
of $C_{\pm}^{i}$ inside $\mathbb{T}^{2}$, with $i=1,\dots,m$, such that system
$Y_{\epsilon}$ produces exactly $1$ limit cycle in $\mathcal{U}_{\pm}^{i}$,
denoted by $\mathcal{O}_{\pm}^{i}$, for each $\epsilon\in]0,\epsilon_{0}]$.
Each limit cycle $\mathcal{O}_{\pm}^{i}$ is hyperbolically attracting and non-
canard ($-$) or repelling and canard ($+$), turns around the unknotted torus
$\mathbb{T}^{2}$ $k$ times vertically and $l$ times horizontally, and tends in
Hausdorff sense to $C_{\pm}^{i}$ in the limit $\epsilon\to 0$. Moreover,
fixing $\epsilon\in]0,\epsilon_{0}]$, for any $\tau\in\mathbb{T}^{2}$ such
that $\tau\notin\mathcal{O}_{\pm}^{i}$, we have that $\alpha(\tau)$ is one of
the repelling cycles $\mathcal{O}_{+}^{i}$ and $\omega(\tau)$ is one of the
attracting cycles $\mathcal{O}_{-}^{i}$.

Theorem 1.1 will follow from Theorem 2.2 of Section 2.2 stated in a more
general framework.

######  Remark 1.

  * •

For $k,l>0$ relatively prime, notice that the slow flow of $Y_{\epsilon}$,
that is $\lim_{\epsilon\to
0}\frac{1}{\epsilon}Y_{\epsilon}|_{\mathcal{O}^{i}_{\pm}}$, is a translation
on $\mathbb{T}^{2}$. It then follows from Theorem 1.1 (and e.g. [2]) that the
limit cycles of $Y_{\epsilon}$ have rotation number $\frac{l}{k}$ and that
$Y_{\epsilon}$ has no dense orbits.

  * •

Related to the previous point, let $J$ denote the horizontal interval
$J=[0,2\pi[\times\left\\{0\right\\}$ and let $\Pi:J\to J$ be the Poincaré map
induced by the flow of $Y_{\epsilon}$ for $\epsilon>0$ sufficiently small. It
follows that $\Pi$ is well-defined, particularly since $\dot{y}=\epsilon$.
This observation, in principle, would allow us to instead consider
$\epsilon$-perturbations of translations over $\frac{l}{k}$ on
$\mathbb{T}^{1}=\mathbb{S}^{1}$. We prefer to not consider this route because
in our general result of Theorem 2.2 this Poincaré map is not necessarily
well-defined as we do not restrict the sign of the slow flow. For example, for
a $(k,l)$-torus knot, a Poincaré map on the torus would need to be defined as
the first return map _after $k$ vertical rotations_. So, if the flow on the
critical manifolds have distinct directions, there would be points on $J$ that
return to $J$ before turning vertically along the torus.

A simple closed curve (i.e. an embedding $\mathbb{S}^{1}\to\mathbb{T}^{2}$)
that turns around the torus $k$ times vertically (meridionally) and $l$ times
horizontally (longitudinally) is often called a $(k,l)$-torus knot. We call a
disjoint finite union of such torus knots a torus link (see Appendix A and [1,
6, 8]). Using this terminology and Theorem 1.1, we can say that for each small
$\epsilon>0$ system $Y_{\epsilon}$ has a link consisting of $2m$ limit cycles
(each limit cycle is a $(k,l)$-torus knot). When $k=1$, $l=0$ or $l=1$, we
deal with trivial knots or unknotted circles (see, for example, system
$X_{\epsilon}$). If $(k,l)=(3,2)$ or $(k,l)=(2,3)$, then the limit cycles of
$Y_{\epsilon}$ are the trefoil knots (see Figure 3). It is not possible to
untangle the trefoil knot into the unknotted circle through
$\mathbb{R}^{3}$-space without “cutting” or “gluing”. If $(k,l)=(5,2)$ or
$(k,l)=(2,5)$, then $2m$ Solomon’s seal knots occur in $Y_{\epsilon}$ for each
$\epsilon>0$ (see Figure 3). For more torus knots see a list in [8].

Figure 3. Examples of knotted critical manifolds of (2) on $\mathbb{T}^{2}$
(left) and on $[0,2\pi]^{2}$ (right) for $m=1$: first row: $(3,2)$-knot
(trefoil); second row: $(2,3)$-knot; third row $(5,2)$-knot (Solomon’s seal);
and fourth row: $(2,5)$-knot. Notice that $(k,l)$-knots are ambient isotopic
(see Appendix A) to $(l,k)$-knots.

In [4] the cyclicity of planar common slow-fast cycles has been studied. The
common slow-fast cycles contain either attracting or repelling portions of the
critical manifold. The authors were focused on the local study of a single
common slow-fast cycle in the plane. In our paper we focus on the global study
of slow-fast vector fields defined on $\mathbb{T}^{2}$ in the presence of a
disjoint collection of common limit periodic sets of torus-knot type.

In Section 2 we define our slow-fast model on the $2$-torus and state the main
result (Theorem 2.2). Section 3 is devoted to the proof of Theorem 2.2. In
Section 4 we explain how we obtain a similar result to Theorem 2.2 in the
presence of regular nilpotent contact points. In Section A we give some basic
definitions and results about torus knots.

##  2\. Assumptions and statement of results

###  2.1. Definition of a slow-fast model on $\mathbb{T}^{2}$ and assumptions

Suppose that $X_{\epsilon,\rho}:\mathbb{T}^{2}\to T\mathbb{T}^{2}$ is a
($C^{\infty}$) smooth $(\epsilon,\rho)$-family of vector fields defined on the
torus $\mathbb{T}^{2}$ of class $C^{\infty}$ where $\epsilon\geq 0$ is the
singular parameter and $\rho\sim\rho_{0}\in\mathbb{R}^{p}$. The tangent bundle
of $\mathbb{T}^{2}$ is denoted by $T\mathbb{T}^{2}$. The parameter $\rho$ is
included for the sake of generality. The first assumption deals with the
dynamics of the fast subsystem $X_{0,\rho}$.

######  Assumption 1.

The system $X_{0,\rho}$ has a smooth $\rho$-family of $2m$ disjoint embedded
closed curves $C_{\rho,-}^{1}$, $C_{\rho,+}^{1}$, $\dots$, $C_{\rho,-}^{m}$,
$C_{\rho,+}^{m}$ of singularities of $X_{0,\rho}$, for some $m\in\mathbb{N}$.
Each curve wraps $k$ times vertically around the torus and $l$ times
horizontally before it comes back to its initial point, for some fixed
relatively prime non-negative integers $k$ and $l$ with $k+l>0$. Moreover,
$C_{\rho,-}^{i}$ (resp. $C_{\rho,+}^{i}$) consists of normally attracting
(resp. repelling) singularities for each $i=1,\dots,m$.

We call the critical manifold of $X_{0,\rho}$ (i.e., the disjoint collection
of $2m$ $(k,l)$-torus knots from Assumption 1) the $2m$-link and denote it by
$C_{\rho}$.

######  Remark 2.

Notice that for $m\geq 2$, the vector field $X_{0,\rho}$ is well-defined only
if the critical manifolds $C_{\rho,\pm}^{i}$, $i=1,\ldots,m$, of $X_{0,\rho}$
alternate each other stability-wise along the fast fibers.

######  Remark 3.

It is worth noting that torus knots and slow-fast dynamics interact in a non-
trivial way. To be more precise, while the two pairs of torus knots of Figure
4 are equivalent (up to homeomorphism and even ambient isotopy) on
$\mathbb{T}^{2}$ [8], they lead to nonequivalent slow-fast systems, as in
Figure 5. In our main Theorem 2.2 we restrict to normally hyperbolic critical
manifolds. This has the advantage of making the proof more concise. However,
as we argue in Section 4, it is possible to extend the results of Theorem 2.2
to some cases where the critical manifold is not normally hyperbolic.

$x$$y$$x$$y$$C_{1}$$C_{2}$$C_{1}^{\prime}$$C_{2}^{\prime}$ Figure 4. An
example of a pair of equivalent (up to homeomorphism and even ambient isotopy)
knots, that is, $C_{i}\sim C_{i}^{\prime}$, $i=1,2$. Compare with Figure 5
$x$$y$$x$$y$$C_{1}$$C_{2}$$C_{1}^{\prime}$$C_{2}^{\prime}$ Figure 5. Even if
the critical manifolds $C_{i}$ on the left and $C_{i}^{\prime}$ on the right
are equivalent as torus knots (compare with Figure 4), they lead to completely
different, and in fact nonequivalent, slow-fast dynamics. In the main part of
this paper we assume that every knot of the $2m$-link critical manifold is
normally hyperbolic, but refer also to Section 4.

A simple topological argument on $\mathbb{T}^{2}$ implies that disjoint torus
knots have to be of the same type (see [8]). More precisely, we have

######  Lemma 2.1.

If $C$ and $D$ are two disjoint knots on $\mathbb{T}^{2}$ of types
$(k_{1},l_{1})$ and $(k_{2},l_{2})$ where $(k_{1},l_{1})$ and $(k_{2},l_{2})$
have the property given in Assumption 1, then $k_{1}=k_{2}$ and $l_{1}=l_{2}$.

In order to ensure that the slow divergence integral of $X_{\epsilon,\rho}$ is
finite (i.e. well-defined), we have to assume that the slow dynamics of
$X_{\epsilon,\rho}$ along the $2m$-link $C_{\rho}$ is regular. Let us recall
that the slow dynamics or the slow flow, denoted by $\widetilde{X}_{\rho}$, is
a $C^{\infty}$-smooth $\rho$-family of one-dimensional vector fields that
describes the passage near $C_{\rho,-}^{i}$ or $C_{\rho,+}^{i}$, with
$i=1,\dots,m$, when $\epsilon$ is positive and small. We have

$\widetilde{X}_{\rho}(\tau)=\lim_{\epsilon\to
0}\frac{(X_{\epsilon,\rho}+0\frac{\partial}{\partial\epsilon})(\tau)|_{M_{\tau}}}{\epsilon}$

where $M_{\tau}$ is any $C^{n}$ center manifold of
$X_{\epsilon,\rho}+0\frac{\partial}{\partial\epsilon}$ at the point $\tau\in
C_{\rho}$ ($n\in\mathbb{N}$ can be as large as we want). For more details
about the definition of the slow dynamics see e.g. [4].

######  Assumption 2.

The slow dynamics $\widetilde{X}_{\rho_{0}}$ is nonzero at each point $\tau\in
C_{\rho_{0}}$.

######  Remark 4.

From Assumption 2 follows that along each of the $2m$ $(k,l)$-torus knots
defined in Assumption 1 the slow dynamics $\widetilde{X}_{\rho_{0}}$ is
regular (thus, without singularities). The slow dynamics may have different
directions along the torus knots. For example, if we replace
$\dot{y}=\epsilon$ in (1) with $\dot{y}=\epsilon\cos(y-x)$, then the slow flow
along $C_{-}$ points upwards and along $C_{+}$ downwards.

Now we can define the slow divergence integral of $X_{\epsilon,\rho}$ along
$C_{\rho,-}^{i}$ (resp. $C_{\rho,+}^{i}$) as:

$I_{-}^{i}(\rho)=\int_{C_{\rho,-}^{i}}\operatorname{div}X_{0,\rho}ds\
\Big{(}\text{resp.
}I_{+}^{i}(\rho)=\int_{C_{\rho,+}^{i}}\operatorname{div}X_{0,\rho}ds\Big{)}$

where $i=1,\dots,m$ and $s$ is the slow time of $\widetilde{X}_{\rho}$. The
slow divergence integral $I_{\pm}^{i}$ is independent of the local chart and
the chosen volume form on $\mathbb{T}^{2}$, and represents the leading order
term of the integral of divergence along orbits of $X_{\epsilon,\rho}$,
multiplied by $\epsilon$ (see Section 3). We typically compute $I_{-}^{i}$ by
dividing the compact curve $C_{\rho,-}^{i}$ into a finite number of segments
$[\tau_{1},\tau_{2}],[\tau_{2},\tau_{3}],\dots,[\tau_{r-1},\tau_{r}]$, where
$r\in\mathbb{N}$, $\tau_{1},\dots,\tau_{r}\in C_{\rho,-}^{i}$ and
$\tau_{1}=\tau_{r}$ ($C_{\rho,-}^{i}$ is closed for all $\rho\sim\rho_{0}$),
and by calculating the integral on each segment in suitable normal form
coordinates:

(3)
$I_{-}^{i}(\rho)=\int_{\tau_{1}}^{\tau_{2}}\operatorname{div}X_{0,\rho}ds+\int_{\tau_{2}}^{\tau_{3}}\operatorname{div}X_{0,\rho}ds+\cdots+\int_{\tau_{r-1}}^{\tau_{r}}\operatorname{div}X_{0,\rho}ds.$

We assume that the points $\tau_{1},\tau_{2},\dots$ follow the direction of
the slow flow along $C_{\rho,-}^{i}$ and that the above normally attracting
segments are small enough such that on each segment we can use a Takens normal
form for $C^{n}$-equivalence

(4) $\left\\{\begin{array}[]{rcl}\dot{x}&=&-x\\\
\dot{y}&=&\epsilon\end{array}\right.$

(see [5]). The slow segment is lying inside $\\{x=0\\}$. Since the divergence
of the vector field in (4) for $\epsilon=0$ is $-1$, it is clear that each
integral in (3) is negative. Thus, the slow divergence integral $I_{-}^{i}$ is
negative for all $\rho\sim\rho_{0}$ and $i=1,\dots,m$. Similarly, we see that
$I_{+}^{i}$ is positive, for all $\rho\sim\rho_{0}$ and $i=1,\dots,m$, because
each $C_{\rho,+}^{i}$ is repelling.

###  2.2. Statement of results

In this section we state our main result. Let $X_{\epsilon,\rho}$ satisfy
Assumptions 1–2. Then for all $\epsilon\sim 0$ and $\epsilon>0$ a $2m$-link of
type $(k,l)$ occurs inside $X_{\epsilon,\rho}$. More precisely,

######  Theorem 2.2.

Suppose that system $X_{\epsilon,\rho}$ satisfies Assumptions 1–2. Then there
exist $\epsilon_{0}>0$, a neighborhood $\mathcal{V}$ of $\rho_{0}$ and tubular
neighborhoods $\mathcal{U}_{\pm}^{i}$ of $C_{\rho_{0},\pm}^{i}$, for each
$i=1,\dots,m$, such that $X_{\epsilon,\rho}$ has exactly one limit cycle in
$\mathcal{U}_{\pm}^{i}$, denoted by $\mathcal{O}_{\epsilon,\rho,\pm}^{i}$, for
all $(\epsilon,\rho)\in]0,\epsilon_{0}]\times\mathcal{V}$ and $i=1,\dots,m$.
The limit cycle $\mathcal{O}_{\epsilon,\rho,-}^{i}$ (resp.
$\mathcal{O}_{\epsilon,\rho,+}^{i}$) is a hyperbolic and attracting non-canard
(resp. repelling canard) limit cycle. Moreover, each
$\mathcal{O}_{\epsilon,\rho,\pm}^{i}$ tends in Hausdorff sense to the
$(k,l)$-torus knot $C_{\rho_{0},\pm}^{i}$ as $(\epsilon,\rho)\to(0,\rho_{0})$.
If $\tau$ is any point on $\mathbb{T}^{2}$ not lying in
$\mathcal{O}_{\epsilon,\rho,\pm}^{i}$, then the $\omega$-limit (resp. the
$\alpha$-limit) of $\tau$ is one of the attracting (resp. repelling) limit
cycles
$\mathcal{O}_{\epsilon,\rho,-}^{1},\dots,\mathcal{O}_{\epsilon,\rho,-}^{m}$
(resp.
$\mathcal{O}_{\epsilon,\rho,+}^{1},\dots,\mathcal{O}_{\epsilon,\rho,+}^{m}$).

We prove Theorem 2.2 in Section 3.

##  3\. Proof of Theorem 2.2

This section is devoted to the proof of Theorem 2.2. Let $X_{\epsilon,\rho}$
satisfy Assumptions 1–2. We divide the proof of Theorem 2.2 into three parts:

$(a)$$\tau_{j}$$\tau_{j+1}$$\tau_{j}$$\tau_{j+1}$$(b)$

Figure 6. A flow box neighborhood with the inset (red) and the outset (blue).
(a) $\epsilon=0$ (b) $\epsilon>0$.

  1. (1)

(Flow-box neighborhoods near $C_{\rho_{0},\pm}^{i}$) Let
$i\in\\{1,\dots,m\\}$. We construct a succession of flow box neighborhoods
along the attracting closed curve $C_{\rho_{0},-}^{i}$ (we can do the same
with the repelling slow curve $C_{\rho_{0},+}^{i}$ by reversing time). Since
$X_{\epsilon,\rho}$ has no singularities on $\mathbb{T}^{2}$ for all small
$\epsilon>0$ and $\rho\sim\rho_{0}$ (Assumptions 1–2 are satisfied), we can
cover $C_{\rho_{0},-}^{i}$ with a finite number of flow boxes that are uniform
in $\epsilon>0$ and $\rho\sim\rho_{0}$, i.e. their size is fixed as
$(\epsilon,\rho)\to(0,\rho_{0})$. (Such flow box neighborhoods don’t exist in
the limit $\epsilon=0$.) This together with a fixed point theorem for segments
will enable us to show the existence of a closed orbit in an
$(\epsilon,\rho)$-uniform tubular neighborhood $\mathcal{U}_{-}^{i}$ of
$C_{\rho_{0},-}^{i}$, for each small $\epsilon>0$ and $\rho\sim\rho_{0}$. Our
approach is based on the techniques from [4].

  2. (2)

(The slow divergence integral along $C_{\rho_{0},\pm}^{i}$) We relate the slow
divergence integral defined along $C_{\rho_{0},-}^{i}$ (Section 2) to the
integral of divergence along orbits of $X_{\epsilon,\rho}$ inside
$\mathcal{U}_{-}^{i}$. Using a well-known connection between the derivative of
the Poincaré map and the divergence integral, and the fact that the slow
divergence integral is negative, we show that $X_{\epsilon,\rho}$ has at most
one limit cycle in $\mathcal{U}_{-}^{i}$ up to shrinking $\mathcal{U}_{-}^{i}$
if needed. The limit cycle is a hyperbolically attracting $(k,l)$-torus knot.
We also use a generalization of Poincaré-Bendixon Theorem to $\mathbb{T}^{2}$
due to Schwartz [13].

  3. (3)

(Global dynamics of $X_{\epsilon,\rho}$ on $\mathbb{T}^{2}$) To study the
global dynamics on $\mathbb{T}^{2}$, we use the same result due to Schwartz
[13].

1\. Flow-box neighborhoods near $C_{\rho_{0},\pm}^{i}$. Let
$\tau\subset\mathbb{T}^{2}$ and let $X$ be a $C^{\infty}$-smooth vector field
on $\mathbb{T}^{2}$. We say that a neighborhood $\mathcal{U}$ of $\tau$ is a
flow box neighborhood of $\tau$ for $X$ if $\mathcal{U}=F([0,1]\times[0,1])$
where $F:[0,1]\times[0,1]\to\mathbb{T}^{2}$ is a $C^{\infty}$-smooth
diffeomorphism with the following property: the vector field $X$ is directed
from outside to inside $\mathcal{U}$ along $F([0,1]\times\\{0\\})$ (the
inset), from inside to outside $\mathcal{U}$ along $F([0,1]\times\\{1\\})$
(the outset), $F(\\{0,1\\}\times[0,1])$ are parts of orbits of $X$ and $X$ has
no singularities in $\mathcal{U}$. Thus, all orbits of $X$ starting at the
inset reach the outset. The set $\tau$ can be a point, a segment, etc. We use
this definition to introduce a slow-fast family of flow box neighborhoods of
$\tau$ for $X_{\epsilon,\rho}$, i.e. a family
$\\{\mathcal{U}_{\epsilon,\rho};(\epsilon,\rho)\in[0,\epsilon_{0}]\times\mathcal{V}\\}$
of neighborhoods of $\tau$ where $\epsilon_{0}>0$, $\mathcal{V}$ is a
neighborhood of $\rho_{0}$,
$\mathcal{U}_{\epsilon,\rho}=F_{\epsilon,\rho}([0,1]\times[0,1])$,
$F_{\epsilon,\rho}:[0,1]\times[0,1]\to\mathbb{T}^{2}$ is an
$(\epsilon,\rho)$-family of smooth diffeomorphisms such that
$\mathcal{U}_{\epsilon,\rho}$ is a flow box neighborhood of $\tau$ for
$X_{\epsilon,\rho}$ for each
$(\epsilon,\rho)\in]0,\epsilon_{0}]\times\mathcal{V}$ and such that the
intersection of the sets $\mathcal{U}_{\epsilon,\rho}$, with
$(\epsilon,\rho)\in[0,\epsilon_{0}]\times\mathcal{V}$, is a neighborhood of
$\tau$. For more details see Definition 9 in [4].

Figure 7. A trefoil knot covered with flow box neighborhoods.

Fix $i\in\\{1,\dots,m\\}$. We divide the compact curve $C_{\rho_{0},-}^{i}$
into the segments
$[\tau_{1},\tau_{2}],[\tau_{2},\tau_{3}],\dots,[\tau_{r-1},\tau_{r}]$ such
that $\tau_{j}\in C_{\rho_{0},-}^{i}$ for all $j=1,\dots,r$,
$\tau_{1}=\tau_{r}$ and such that the points
$\tau_{1},\tau_{2},\dots,\tau_{r}$ follow the direction of the slow flow along
$C_{\rho_{0},-}^{i}$ (i.e. the slow flow of $X_{0,\rho_{0}}$ goes from
$\tau_{1}$ to $\tau_{2}$, from $\tau_{2}$ to $\tau_{3}$, etc.). For each
segment $[\tau_{j},\tau_{j+1}]\subset C_{\rho_{0},-}^{i}$, $j=1,\dots,r-1$, we
define a slow-fast family of flow box neighborhoods of $[\tau_{j},\tau_{j+1}]$
inside $X_{\epsilon,\rho}$ (see Figure 6). We describe it using equivalence
normal form coordinates. There exists a local chart on $\mathbb{T}^{2}$ around
$[\tau_{j},\tau_{j+1}]$ in which $X_{\epsilon,\rho}$ is locally given by the
Takens normal form (4) (the segment $[\tau_{j},\tau_{j+1}]$ can be as small as
we need). In the normal form coordinates, the segment is given by $[-1,1]$ on
the $y$-axis. The outset is a parabola like segment above the line
$\\{y=1\\}$. From the end points of the outset we have two parts of orbits of
$X_{\epsilon,\rho}$ (they are horizontal in the limit $\epsilon=0$). The inset
consists of two lines and a convex segment between them as indicated in Figure
6. Note that the vector field (4) is directed from outside to inside the
neighborhood along the inset for each small $\epsilon>0$ because the slow
dynamics $y^{\prime}=1$ points upwards near the $y$-axis. A detailed
description of the slow-fast family of flow box neighborhoods can be found in
Section 5 of [4].

Thus, we have constructed a finite number of flow box neighborhoods
$\mathcal{U}_{\epsilon,\rho}^{j}$, $j=1,\dots,r-1$, along $C_{\rho_{0},-}^{i}$
(see Figure 7). Since each flow box neighborhood is contracted along its
segment, we may assume that the outset of $\mathcal{U}_{\epsilon,\rho}^{j}$ is
completely inside $\mathcal{U}_{\epsilon,\rho}^{j+1}$, for each
$j=1,\dots,r-2$, and that the outset of $\mathcal{U}_{\epsilon,\rho}^{r-1}$
lies inside $\mathcal{U}_{\epsilon,\rho}^{1}$ (see Proposition 4 of [4]). This
implies that the Poincaré map from the inset of
$\mathcal{U}_{\epsilon,\rho}^{1}$ to itself is well defined and smooth for
$\epsilon>0$ small enough. Since the inset is a segment, there is a fixed
point by Brouwer’s Fixed Point Theorem. Thus, there exists a closed orbit near
$C_{\rho_{0},-}^{i}$ for each small $\epsilon>0$ and $\rho\sim\rho_{0}$.

2\. The slow divergence integral along $C_{\rho_{0},\pm}^{i}$. Since
$C_{\rho_{0},-}^{i}$ is normally attracting and compact, it is well known (see
e.g. [5]) that for any small $\kappa>0$ we can find $\epsilon_{0}>0$, a
neighborhood $\mathcal{V}$ of $\rho_{0}$ and a tubular neighborhood
$\mathcal{U}_{-}^{i}$ of $C_{\rho_{0},-}^{i}$ such that

(5)
$\frac{I_{-}^{i}(\rho_{0})-\kappa}{\epsilon}\leq\int_{\mathcal{O}_{\epsilon,\rho}}\operatorname{div}{X_{\epsilon,\rho}}dt\leq\frac{I_{-}^{i}(\rho_{0})+\kappa}{\epsilon}$

where $I_{-}^{i}(\rho_{0})<0$ is the slow divergence integral defined in
Section 2 and $\mathcal{O}_{\epsilon,\rho}$ is an arbitrary closed orbit of
$X_{\epsilon,\rho}$ in $\mathcal{U}_{-}^{i}$ with
$\epsilon\in]0,\epsilon_{0}]$ and $\rho\in\mathcal{V}$. We use now (5) to show
that any closed orbit $\mathcal{O}_{\epsilon,\rho}$ of $X_{\epsilon,\rho}$ in
$\mathcal{U}_{-}^{i}$ is hyperbolically attracting. It suffices to take a
$\kappa>0$ such that $\kappa<-I_{-}^{i}(\rho_{0})$ and to use the Poincaré
formula
$P^{\prime}(s_{0})=e^{\int_{\mathcal{O}_{\epsilon,\rho}}\operatorname{div}{X_{\epsilon,\rho}}dt}$
where $P$ is the Poincaré map defined on a transverse section near
$\mathcal{O}_{\epsilon,\rho}$ parametrized by a regular parameter $s$
($\mathcal{O}_{\epsilon,\rho}$ corresponds to $s_{0}$). Note that the torus is
orientable. Since the integral of divergence in (5) is negative (and of order
$O(\frac{1}{\epsilon})$), we have $|P^{\prime}(s_{0})|<1$ and the closed orbit
$\mathcal{O}_{\epsilon,\rho}$ is a hyperbolically attracting limit cycle.

To prove that $C_{\rho_{0},-}^{i}$ can produce at most one limit cycle inside
$\mathcal{U}_{-}^{i}$, we use the following result due to Schwartz [13] (see
also Theorem 6.6 in [3]):

######  Theorem 3.1.

Assume that $x^{\prime}=X(x)$ is an autonomous system of class $C^{2}$ on a
compact, connected two-dimensional orientable manifold $M$ of class $C^{2}$.
Assume that this system defines a flow and that the manifold $M$ is not a
minimal set. Then, if the $\omega$-limit set $\omega(\tau)$ of a point
$\tau\in M$ does not contain any critical point, $\omega(\tau)$ must be
homeomorphic to the circle $\mathbb{S}^{1}$ (i.e. is a periodic trajectory).

We apply Theorem 3.1 to $M=\mathbb{T}^{2}$ and $X=-X_{\epsilon,\rho}$ for each
fixed $\epsilon>0$ and $\rho\sim\rho_{0}$ (in this paper, both
$\mathbb{T}^{2}$ and $X_{\epsilon,\rho}$ are of class $C^{\infty}$). Let us
recall that a set $S\subset M$ is minimal if it is nonempty, invariant, closed
and there is no proper subset of $S$ that has all these properties. Since
$X_{\epsilon,\rho}$ has at least one closed (or periodic) orbit for each
$\epsilon>0$ small enough (see Step 1), $\mathbb{T}^{2}$ is not minimal. As
already mentioned, $X_{\epsilon,\rho}$ has no critical points on
$\mathbb{T}^{2}$ for $\rho\sim\rho_{0}$ and $\epsilon>0$ small enough.

Suppose now that $X_{\epsilon,\rho}$ has at least two limit cycles in
$\mathcal{U}_{-}^{i}$ ($\mathcal{O}_{1}$ and $\mathcal{O}_{2}$). Then
$\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ are $(k,l)$-torus knots and bound an
invariant set $S$. Take any $\tau\in S$ not lying on a closed orbit of
$X_{\epsilon,\rho}$ (such $\tau$ exists if we are close enough to the isolated
closed orbits $\mathcal{O}_{1}$ or $\mathcal{O}_{2}$). Following Theorem 3.1,
$\omega(\tau)\subset S$ w.r.t. $-X_{\epsilon,\rho}$ is a periodic trajectory
different from $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ ($\mathcal{O}_{1}$ and
$\mathcal{O}_{2}$ are repelling for $-X_{\epsilon,\rho}$). This leads to a
contradiction because the periodic trajectory $\omega(\tau)$ is not repelling
(all closed orbits generated by $C_{\rho_{0},-}^{i}$ are hyperbolic and
repelling limit cycles w.r.t. $-X_{\epsilon,\rho}$).

From Step 1 and Step 2 follows that there exists $\epsilon_{0}>0$, a small
neighborhood $\mathcal{V}$ of $\rho_{0}$ and a tubular neighborhood
$\mathcal{U}_{-}^{i}$ of $C_{\rho_{0},-}^{i}$ such that for each
$(\epsilon,\rho)\in]0,\epsilon_{0}]\times\mathcal{V}$ system
$X_{\epsilon,\rho}$ has precisely one limit cycle in $\mathcal{U}_{-}^{i}$. It
is hyperbolic and attracting and clearly of $(k,l)$-torus knot type. By
reversing time, we can prove a similar result near $C_{\rho_{0},+}^{i}$ for
each $i=1,\dots,m$ ($C_{\rho_{0},+}^{i}$ generates one hyperbolic and
repelling limit cycle).

3\. Global dynamics of $X_{\epsilon,\rho}$ on $\mathbb{T}^{2}$. First note
that the $2m$ limit cycles obtained in Step 1 and Step 2 are the only possible
periodic trajectories of $X_{\epsilon,\rho}$. Indeed, any orbit with the
initial point
$\tau\in\mathbb{T}^{2}\setminus\cup_{i=1}^{m}\mathcal{U}_{\pm}^{i}$ ($\tau$ is
thus uniformly away from the critical manifold of $X_{0,\rho_{0}}$) cannot be
periodic (it is attracted to an attracting closed part $C_{\rho_{0},-}^{i}$
and stays there forever). Theorem 3.1 implies now the rest of Theorem 2.2.

##  4\. Regular nilpotent contact points of finite order

In this section we briefly argue that a similar result to Theorem 2.2 can be
obtained, without much effort, in some cases where the critical manifold
$C_{\rho}$ is not normally hyperbolic. Instead of providing all the
technicalities, we refer to [4], and only indicate the main ideas.

######  Definition 1.

A point $\tau\in C_{\rho_{0}}$ is called a _nilpotent contact point_ if the
linear part of $X_{0,\rho_{0}}$ at $\tau$ (computed in local coordinates) is
nilpotent. Near such a point, with $(\epsilon,\rho)\sim(0,\rho_{0})$, system
$X_{\epsilon,\rho}$ is $C^{\infty}$-equivalent to
$\\{\dot{x}=y-f(x,\rho),\dot{y}=\epsilon(g(x,\epsilon,\rho)+O(y-f(x,\rho))\\}$
where $(x,y)$ are local normal form coordinates, $\tau$ is given by
$(x,y)=(0,0)$, $f,g,O$ are $C^{\infty}$-smooth and the $1$-jet
$j_{1}f(x,\rho_{0})$ at $x=0$ is zero (see [4]). A nilpotent contact point
$\tau$ is further called _regular_ if $g(0,0,\rho_{0})\neq 0$. Finally, a
nilpotent contact point $\tau$ is called of finite order, if there exists an
integer $n\geq 2$ such that the vanishing order of $f(x,\rho_{0})$ at $x=0$ is
equal to $n$.

It is not difficult to see that the above definition is independent of the
chosen normal form.

It turns out that the main ideas of Theorem 2.2 also hold to the case where
the critical manifold may contain isolated regular nilpotent contact points of
finite order. To be more precise, Assumption 1 can be relaxed in the following
sense:

Assumption 1’:

The system $X_{0,\rho}$ has a smooth $\rho$-family of $2m$ disjoint embedded
curves $C_{\rho,-}^{1}$, $C_{\rho,+}^{1}$, $\dots$, $C_{\rho,-}^{m}$,
$C_{\rho,+}^{m}$ of singularities of $X_{0,\rho}$, for some $m\in\mathbb{N}$.
Each curve wraps $k$ times vertically around the torus and $l$ times
horizontally before it comes back to its initial point, for some fixed
relatively prime non-negative integers $k$ and $l$ with $k+l>0$. Moreover,
$C_{\rho,\pm}^{i}$ consists of normally hyperbolic singularities for each
$i=1,\dots,m$ except at a finite number of points where $C_{\rho,-\pm}^{i}$
(resp. $C_{\rho,+}^{i}$) may have isolated regular nilpotent contact points of
finite order.

Then, one can show that under Assumptions 1’ and 2 (Assumption 2 is valid at
normally hyperbolic singularities and the slow flow is unbounded at the
contact points), $X_{\epsilon,\rho}$ has exactly $2m$ limit cycles from which
$m$ are attracting and $m$ are repelling. The proof would follow the same
arguments as in Section 3 and [4]. Instead of describing these technicalities,
we simply recall that the passage near regular nilpotent contact points of
finite order is a local contraction, see the details in [4] where also Figure
5 is insightful. Examples of knotted critical manifolds with regular nilpotent
contact points of finite order are provided in Figures 8 and 9.

Figure 8.  Example of a $2$-link critical manifold consisting of two
$(1,1)$-knots with an odd regular contact point. An odd contact point is
locally given by the singularity $y=x^{2n+1}$, $n\in\mathbb{N}$. On the left
we sketch the singular limit, observe that the critical manifold has an odd
contact with the fast foliation at $x=\pi$. The red knot corresponds to the
attracting critical manifold, while the blue knot to the repelling one. On the
right we sketch possible orbits for $\epsilon>0$ sufficiently small, where the
red knot indicates the attracting non-canard limit cycle and the blue knot the
repelling canard limit cycle. Orbits away from the canard limit cycle approach
the non-canard limit cycle as $t\to\infty$. Figure 9.  Example of a $2$-link
critical manifold consisting of two Trefoils ($(3,2)$-knots) with regular
contact points. On the left we sketch the singular limit, observe that the
critical manifolds (thick red and blue knots) contain, for example, jump
points. All flows are shown in forward time. The shaded red path indicates a
singular slow-fast cycle, that is a knot consisting of slow and fast orbits of
the singular limit. As sketched on the right, the singular slow-fast knot
perturbs to a non-canard attracting limit cycle, which in fact corresponds to
a relaxation oscillation. On the other hand, the blue knot perturbs to a
repelling canard limit cycle, also of slow-fast type (singular limit not
shown). The main argument to show that a unique limit cycle bifurcates from
each singular slow-fast knot for $\epsilon>0$ sufficiently small, is that the
local passage through the regular contact points is a contraction.

Although we do not provide a specific model of a smooth slow-fast system on
$\mathbb{T}^{2}$ with knotted critical manifold with contact points, we point-
out that it is possible to construct such systems using, for example, Hermite
interpolation [14].

##  Appendix A Basic notions of Knot Theory

In this section we present some basic notions, definitions, and results about
knot theory and, in particular, torus knots. For more details see, for
example, [1, 8]. We start with the definition of a knot:

######  Definition 2.

A _knot_ is an embedding $K:\mathbb{S}^{1}\to\mathbb{R}^{3}$ of a $1$-sphere
into $\mathbb{R}^{3}$. A _link_ $L:\cup\,\mathbb{S}^{1}\to\mathbb{R}^{3}$ is a
disjoint, finite union of knots222More generally, one can define a knot as
being an embedding $K:\mathbb{S}^{r}\to X$, where $X$ is a topological space
(usually $\mathbb{S}^{n}$ or $\mathbb{R}^{n}$). Alternatively, one can
consider a knot as being a subset homemorphic to $\mathbb{S}^{r}$, $r\geq 1$..

In order to introduce a notion of equivalence between knots, we have the
following:

######  Definition 3.

A homotopy $h_{t}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ is called an _ambient
isotopy_ if $h_{0}$ is the identity and every $h_{t}$ is a homeomorphism.

There are several notions of equivalence for knots. For our purposes it will
suffice to keep in mind the following two:

######  Definition 4.

  1. (1)

Two knots (or links) $K$ and $K^{\prime}$ are _equivalent_ , if there is a
homeomorphism $h:\mathbb{R}^{3}\to\mathbb{R}^{3}$ such that $h(K)=K^{\prime}$.

  2. (2)

Two knots (or links) $K$ and $K^{\prime}$ are _ambient isotopic_ , if there is
an ambient isotopy such that $h_{1}(K)=K^{\prime}$.

In this paper we are concerned with torus knots, that is, knots that lie in
the $2$-torus
$\mathbb{T}^{2}=\mathbb{S}^{1}\times\mathbb{S}^{1}\subset\mathbb{R}^{3}$. In
this regard, we say that a torus knot $K$ is of class $(k,l)$ if the knot
winds the torus $k$-times vertically and $l$-times horizontally, and where $k$
and $l$ are coprime (see e.g. Figure 3). It turns out that the only
nonequivalent knots (up to homeomorphism) are those presented in Figure 10.

$K_{1}$$K_{2}$ Figure 10. Representatives of the two classes of nonequivalent
torus knots (up to homeomorphism). That is, every knot on the $2$-torus can be
transformed via a homeomorphism to either $K_{1}$, which is a knot of class
$(1,0)$, or $K_{2}$, which is a knot of class $(0,0)$ [8].

We see that equivalence up to homeomorphism is too coarse. In contrast,
ambient isotopies preserve orientation, which leads to a finer classification.

######  Theorem A.1 ([6, 8]).

Let $K_{i}$, $i=1,2$, be torus knots of class $(k_{i},l_{i})$ respectively.
Then $K_{1}$ and $K_{2}$ are ambient isotopic if and only if
$(k_{1},l_{1})=\pm(k_{2},l_{2})$. If one considers only the case where $k_{i}$
and $l_{i}$ are positive, then $K_{1}$ and $K_{2}$ are ambient isotopic if and
only if $k_{1}=k_{2}$ and $l_{1}=l_{2}$ or $k_{1}=l_{2}$ and $l_{1}=k_{2}$.

Examples of equivalent knots modulo ambient isotopy can be seen in Figure 3.

## References

  * [1] C. C. Adams.  The knot book.  American Mathematical Society, Providence, RI, 2004.  An elementary introduction to the mathematical theory of knots, Revised reprint of the 1994 original.
  * [2] H. Broer and F. Takens.  Dynamical systems and chaos, volume 172.  Springer Science & Business Media, 2010.
  * [3] K. Ciesielski.  The Poincaré-Bendixson theorem: from Poincaré to the XXIst century.  Cent. Eur. J. Math., 10(6):2110–2128, 2012.
  * [4] P. De Maesschalck, F. Dumortier, and R. Roussarie.  Cyclicity of common slow-fast cycles.  Indag. Math. (N.S.), 22(3-4):165–206, 2011.
  * [5] F. Dumortier and R. Roussarie.  Canard cycles and center manifolds.  Mem. Amer. Math. Soc., 121(577):x+100, 1996.  With an appendix by Li Chengzhi.
  * [6] R. W. Ghrist, P. J. Holmes, and M. C. Sullivan.  Knots and links in three-dimensional flows, volume 1654 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 1997.
  * [7] J. Guckenheimer and Yu. Ilyashenko.  The duck and the devil: canards on the staircase.  Mosc. Math. J., 1(1):27–47, 2001.
  * [8] D. Rolfsen.  Knots and links.  Publish or Perish, Inc., Berkeley, Calif., 1976.  Mathematics Lecture Series, No. 7.
  * [9] I. Schurov and N. Solodovnikov.  Duck factory on the two-torus: multiple canard cycles without geometric constraints.  J. Dyn. Control Syst., 23(3):481–498, 2017.
  * [10] I. V. Schurov.  Canard cycles in generic fast-slow systems on the torus.  Tr. Mosk. Mat. Obs., 71:200–234, 2010.
  * [11] I. V. Schurov.  Ducks on the torus: existence and uniqueness.  J. Dyn. Control Syst., 16(2):267–300, 2010.
  * [12] I.V. Schurov.  Duck farming on the two-torus: multiple canard cycles in generic slow-fast systems.  Discrete Contin. Dyn. Syst., (Dynamical systems, differential equations and applications. 8th AIMS Conference. Suppl. Vol. II):1289–1298, 2011\.
  * [13] A. J. Schwartz.  A generalization of a Poincaré-Bendixson theorem to closed two-dimensional manifolds.  Amer. J. Math. 85 (1963), 453-458; errata, ibid, 85:753, 1963.
  * [14] A Spitzbart.  A generalization of Hermite’s interpolation formula.  The American Mathematical Monthly, 67(1):42–46, 1960.


# Thermo-optical reshaping of SHG emission from dimer all-dielectric
nanoresonators

Olesia Pashina  Physics and Engineering Department, ITMO University, St.
Petersburg, Russia, 197101  Kristina Frizyuk  Physics and Engineering
Department, ITMO University, St. Petersburg, Russia, 197101  George Zograf
Physics and Engineering Department, ITMO University, St. Petersburg, Russia,
197101  Mihail Petrov  Physics and Engineering Department, ITMO University,
St. Petersburg, Russia, 197101  Corresponding author<EMAIL_ADDRESS>
###### Abstract

In recent years resonant semiconductor and all-dielectric nanophotonics
offered a lot of possibilities for thermally-induced light manipulation at the
nanoscale. Owing to high-quality resonant states, such nanostructures allow
for efficient light-to-heat conversion supported by various temperature
detection approaches based on thermally sensitive intrinsic optical responses
(photoluminescence, Raman scattering, thermorefraction etc.). In this work, we
study theoretically a phenomenon of the photothermal reshaping of the
radiation pattern of second-harmonic generation (SHG) that occurs in resonant
all-dielectric systems. In the suggested geometry near-IR pulsed laser is
utilized for excitation of SHG while simultaneous continuous wave visible
laser heats the structure. The thermooptical switching of the resonant optical
states in the nanostructures governs the reconfiguration of the emission
pattern, without significant loss in the magnitude of the SHG. We believe,
that our findings will pave the way for subwavelength-size near-IR thermally
switchable nonlinear optical devices.

††journal: ol

After the first observation of the phenomenon of optical harmonics generation
[1], such nonlinear optical effect proved its high efficiency in crystals and
waveguides [2], becoming a commercially available technology for optical
microscopy [3, 4], or laser sources and even green-light laser pointers. While
high efficiencies of SHG in noncentrosymmetric bulk materials could be
achieved by a phase-matching [5], going down to subwavelength-scale particles
and retaining high efficiency was quite a challenge. Different plasmonic [6,
7, 8, 9, 10, 11, 12, 13], all-dielectric [14, 15, 16, 17, 18, 19, 20, 21, 22],
and even hybrid (metal-dielectric) [23, 24, 25, 26, 27, 28, 29] designs were
proposed, however, the conversion of the pump to second harmonic (SH) signal
was still of the order of $10^{-5}-10^{-4}$. Moreover, very recent attractive
family of materials for nanophotonics community - transition metal
dichalcogenides (TMDCs) - both monolayers [30] and nanoresonators [31]
demonstrated its potential in nonlinear nanophotonics. Even though TMDC
monolayers mostly posses a very high nonlinear susceptibility, the bulk
second-order nonlinearity tensor is zero due to crystalline symmetry for most
of the materials [32]. There are other approaches for boosting SH conversion
based on bound states in the continuum (BIC) or quasi-BIC modes in resonant
structures (arrays or metasurfaces) reaching $10^{-2}-10^{-1}$, but these
structures are no longer comparable to the wavelength scale [33, 34, 35].
Moreover, high-Q states in single cylindrical subwavelength resonators
supported $10^{-1}$ SH conversion [36], however tunability of modes and
radiation pattern still remains a challenge.

On the other hand, such nanoantennas proved their ability for strong optical
heating and thermometry due to resonantly enhanced optical absorption [37, 38,
39]. Owing to strong thermooptical effects [40, 41] all-dielectric structures
already demonstrated their proficiency in modulation of linear properties of
propagating [42, 43] and scattered [44] light. Moreover, the possibility of
all-optical thermal modulation of the SHG emission intensity in resonant
semiconductor nanostructures has been reported [45] recently, which makes it a
an accessible experimental technique.

Figure 1: Schematic of the concept. Continuous wave heating of a dimer
nanostructure steers the SHG emission far-field pattern. The dimer consists of
a triangular prism and a cylinder of AlxGa1-xAs located on a glass substrate.

In this work, we propose in situ thermooptical tuning of SHG pattern in bi-
resonant dimer nanoantenna, consisting of two nanoresonators, thereby making
the whole nanoantenna of subwavelength size (see Fig. 1). Being a highly
important property of nanoscale light emission, the directivity pattern of the
SHG was well studied for single resonant nanoantennas [46, 47]. However, under
the intense continuous wave (CW) laser excitation [45, 44] the spectral
position of resonances can shift due to thermooptical effects leading to
reconfiguration of the emission pattern. In this view, we  advance the
approach suggested in [45] proposing two closely located nanostructures
forming a dimer nanoantenna. Different particles shape and slight spectral
detuning between the resonances allow for reconfiguration of the nonlinear
emission pattern with the temperature increase as schematically shown in Fig.
1. Such an approach, when each single element is tuned at particular
temperature regime, can be utilized for desingning thermally stable or tunable
nonlinear optical systems.

Figure 2: Resonant properties of a single nanocylinder for the incident and
second harmonic frequencies. A nanodisk with a diameter
$\mathrm{d}=386\>\mathrm{nm}$ and a height $\mathrm{h}=400\>\mathrm{nm}$,
placed on the glass substrate is excited with normally incident plane wave.
The orange color shows the magnetic dipole resonance of linear scattering, the
violet color refers to the electric octupole resonance of the SHG. The SHG
intensity was normalized to the square of the incident intensity
$\mathrm{I_{SHG}/I_{0}^{2}}$, and the scattering cross-section was normalized
to the geometric section of the particle ${\sigma_{sca}/(\pi R^{2}))}$, where
$R$ is the radius of the cylinder. Far-field radiation patterns are shown for
both linear scattering and SHG. In the lower left corner the polarization
arrows and the magnitude of the Electric field on the surface of the
nanoparticle are shown in green.

We rely on semiconductor solid solution of AlxGa1-xAs which is an efficient
nonlinear material with large quadratic nonlinear optical coefficient
$d\approx 100$ pm/V at wavelengths of about 1550 nm [48, 49].  This
coefficient is related to the second-order nonlinear susceptibility as
$d=\frac{1}{2}\chi^{(2)}$ [32] and determines the nonlinear polarization for
the material under consideration according to the equation:
$P_{i}=2\varepsilon_{0}d\sum_{jk}\delta_{ijk}E_{j}E_{k}$, where $\delta_{ijk}$
vanishes, if any of the two indices i, j, k coincide and $\delta_{ijk}=1$
otherwise. For further simulations, we use particular compound $x=0.198$ which
properties are available elsewhere [50, 51, 52, 53]. We start with analyzing
the nonlinear and thermooptical response of a single cylindrical
nanostructure. The parameters of a cylinder structure can be optimized in a
way to achieve double resonant condition both at the fundamental and SH
wavelengths. Indeed, the resonant mode excitation at the fundamental
wavelength is vital for achieving the enhancement of the SHG signal [54, 29,
55, 56, 57, 58]. On the other hand, the resonance at second harmonic emission
wavelength $\lambda_{\mathrm{pump}}/2$ provides additional enhancement of the
SHG signal and, more importantly, governs the SH emission pattern [46, 54]).
Fig. 2 depicts the elastic light scattering properties and SHG efficiency of a
single AlGaAs nanodisk optimized for pumping in the infrared region 1550 -
1600 nm typical for commercially available pulsed-laser systems. At roughly
1570 nm the nanodisk of 386 nm diameter and 400 nm height possesses resonant
magnetic dipole response, allowing for enhancement of the electromagnetic
field inside the nanodisks volume.  At the same time, at the SH wavelength the
sharp resonant enhancement of the emission is observed. The distribution of
the SH electric field on the nanoparticle surface shown in Fig. 2, as well as
the corresponding radiation pattern with four pronounced lobes due to the
substrate leakage, are characteristic of the octupole mode with the axial
number m=2 [59].  The pronounced lobes in the lower direction correspond to
the preferential substrate emission.

Figure 3: Critical-coupling heating regime and self-consistent optical heating
of a single nanoparticle. a) The ratio of radiation and ohmic losses
$\gamma_{\mathrm{rad}}/\gamma_{\mathrm{Ohm}}$ for the eigenmode of cylinders
with optimized ($\mathrm{h}=380\>\mathrm{nm}$ and
$\mathrm{d}=373\>\mathrm{nm}$) and non-optimized
($\mathrm{h}=400\>\mathrm{nm}$ and $\mathrm{d}=386\>\mathrm{nm}$) size. The
optimized case provides the critical-coupling heating regime, corresponding to
the equality of ohmic and radiation losses
$\gamma_{\mathrm{rad}}=\gamma_{\mathrm{Ohm}}$. b) The temperature increase
$\delta\mathrm{T}$ dependence on the CW laser intensity at optimized (orange
lines) and non-optimized (violet lines) cylinder parameters. The solid lines
show are obtained with simulations without taking into account temperature-
dependent optical constants, while the dashed lines show the self-consistent
optical heating modeling. Figure 4:  Thermally tunable bi-resonant dimer
radiation patterns. a) Schematic representation of a dimer structure
consisting of of a cylinder and a prism with the following geometric
parameters $\mathrm{h}=400\>\mathrm{nm}$, $\mathrm{d}=386\>\mathrm{nm}$,
$\mathrm{a}=350\>\mathrm{nm}$, $\mathrm{b}=460\>\mathrm{nm}$,
$\mathrm{r}=200\>\mathrm{nm}$. The orange semicircle represents the hemisphere
of the second harmonic signal detection from the dimer. b) The dependence of
SHG intensity normalized to the square of the incident intensity
$\mathrm{I_{SHG}/I_{0}^{2}}$, ($\mathrm{I_{0}}=10^{13}\>\mathrm{W/m^{2}}$[60])
at two temperature values for a dimer consisting of a prism and a cylinder.
The blue line corresponds to 25o C and the orange line is for the case of the
elevated temperature of $325^{o}$ C. The fundamental wavelength is 1565 nm and
is marked with the vertical black line. c) Top view of the SHG radiation
detected on the upper hemisphere of a dimer consisting of a cylinder and a
prism at a wavelength of 1565 nm at room $\delta\mathrm{T}=0\>\mathrm{K}$ and
elevated $\delta\mathrm{T}=300\>\mathrm{K}$ temperature. The electric field E
is polarized along the x-axis. The radiation patterns at room temperature and
at elevated temperature have the same colorscale. The dominant directions of
strong switching for SH intensities during heating are marked with white
points 1 and 2, and the detection hemisphere is divided into left and right
parts by a white dotted line.

At the next step, we pay attention to thermooptical response of the considered
structures. In accordance with the recent experimental successes [45], we
consider the normal incidence of continuous wave laser beam that heats the
nanostructure. The laser-assisted optical heating has a big advantage of
tunable and local heating of the nanostructure [61] rather than changing the
temperature of the whole setup. Owing to the strong spectral dependence of the
material losses in the semiconductor structures, the proper choice of the
heating laser wavelength and excitation conditions is crucial on the way of
achieving strong thermooptical response. It is well known that the radiative
($\gamma_{\mathrm{rad}}$) and non-radiative ($\gamma_{\mathrm{Ohm}}$) losses
rate of particular mode should be balanced [37, 62, 38] providing the optimal
heating condition. One can easily balance in semiconductors in the region of
strong interband absorption, which in the case of crystalline silicon occurs
around 700 nm region. Having fixed the ratio of the height and diameter of the
nanocylinder $\mathrm{h/d}$, we slightly scale the geometric parameters of the
resonator. As a result, the eigenmode is spectrally shifted and, while the
radiative losses do not change significantly, the non-raidiative losses are
fastly changing due to their strong spectral dependence (see Fig. 3 (a)).
Therefore, the ratio $\gamma_{\mathrm{rad}}$/$\gamma_{\mathrm{Ohm}}$ takes
different values reaching the optimal condition for disks of 373 nm diameter.
Worth mentioning, that we are neglecting the pulsed laser heating assuming
that the pulse intensity is short enough and the repetition rate is low [63].

The thermooptical simulations performed with commercially available COMSOL
Multyphiscs software show that the heating efficiency of the optimized
structure (orange line) is much higher than that of the non-optimized (blue)
which is a clear sign of critical coupling regime. The intensity range was
chosen well below the damage threshold[60]. However, these data were obtained
for temperature-independent material losses without taking into account the
thermooptical effects. At the same time, intense optical heating induces a
change in the refractive index of the nanoobject material according to a
linear law $n(T)=n_{0}+\alpha T$, where $\alpha$ is the thermorefractive
coefficient. The typical value of thermorefractive coefficients for AlAs
$\alpha=1.25\cdot 10^{-4}{K}^{-1}$, at 750 nm wavelength [51]. The imaginary
part of the refractive index also has strong dependence on temperature $k(T)$.
However, these parameters also depend both on the wavelength and particular
compound of solid mixture $x$, and they were interpolated from the data
available for AlAs and GaAs [64, 65] taken in the appropriate proportion
corresponding to solid mixture AlxGa1-xAs for $x=0.198$. Surprisingly,
including the thermorefractive effect $n(T),k(T)$ into the simulations leads
to rearrangement of the losses channels in the nanoresonator: absorption in
initially non-optimized geometry fastly increases with the temperature and to
match the radiative losses and locking the system in the balance losses
condition. Now, the originally optimized and non-optimized disks have very
similar intensity dependence (Fig. 3b), and non-optimized geometry shows even
slightly stronger thermooptical response, which we utilized for further
simulations. The optimal heating conditions for temperature dependent optical
parameters is a matter of independent study and will be reported in the
following works. The figure depicts a steady-state solution of CW laser
heating of single nanodisk on a glass substrate. We assume the establishment
of average temperature within one nanoparticle and the dimer as whole, since
the thermal conductivity of our semiconductor nanostructure $\kappa_{AlGaAs}$
= 22.75 W m/K is much higher than the thermal conductivity of air
$\kappa_{air}$= 0.026 W m/K and glass $\kappa_{SiO_{2}}$= 0.8 W m/K, adjacent
to the surface of nanoparticles. In this regard, homogeneous thermal
distribution within our nanodimer becomes a fair approximation, because an
extremely small part of the heat leaks into the substrate or into the air [66,
38].

Thus, one of the main thermooptical effects observed in the nanoresonator
structures is related to the spectral shift of the resonant spectral line due
to thermorefractive effect. Inevitably, that leads to taking the system out of
the resonant condition both at fundamental and SH wavelengths. As soon as the
dominant multipole contribution is suppressed the other non-resonant
multipoles start to influence the nonlinear response leading to
reconfiguration of the far-field radiation pattern. In this regard, we
consider a nanoantenna structure compound of two elements (dimer structure),
which resonances are detuned from each other for the value corespondent to the
value of thermal drift of the resonance. Thus, one of the elements of the
dimer will be in the resonant state in the cold regime, while another at the
elevated temperature. Speaking of experimental realization, it was already
experimentally shown that SH emission can be thermally detuned and detected
[45], moreover, asymmetric nanoresonators can be fabricated as well [47]. In
order to achieve reconfiguration and switching of the SH far-field pattern,
one should also choose the different shapes and symmetries of "on" and "off"
structures. In this regard, we have chosen a cylinder-prism structure. One may
fairly ask a question why should it be particularly a prism. The reason for
that lies in the symmetry governed nonlinear response of nanoantenna
structures: by lowering the symmetry of the structure one enables the
excitation of the modes, which are prohibited in the case of a disk [67]. In
this regard, the sphere, cube, or cone have too high symmetries, while prism
has well-pronounced resonances, low symmetry, and are also relatively easy to
fabricate with planar technique. The cylinder-prism structure under the
consideration is shown in Fig. 4 (a) along with the geometry of the excitation
and collection aperture (shown in yellow circle). In the proposed structure,
in accordance with Fig. 5, the cylinder possesses the SH resonance at the
fundamental wavelength and therefore is in the "on"-regime at room
temperature, while prism gets into the resonant state at the elevated
temperature. The simulated SH spectrum from dimer structure is shown in Fig. 4
(b) at room temperature (blue line) and at the elevated temperature of 325o C
(orange line). One can see that the SH spectrum has a two-peak structure,
which are separated by a spectral gap of 50 nm and inherent the spectral
detuning between the resonances of cylinder and prism structures. Once
designed to have the long-wavelength SHG resonance at 1565 nm at room
temperature, the short-wavelength peak gets into the resonant state at
elevated temperature. This is accompanied by strong reshaping of the SH far-
field patterned which is governed by different responses of the cylinder and
prism structures. For the two most dominant directions, we obtained the ratio
of the SH intensities
$\mathrm{\left.I\right|_{325^{o}C}/\left.I\right|_{25^{o}C}}=0.5$ and
$\mathrm{\left.I\right|_{325^{o}C}/\left.I\right|_{25^{o}C}}=23.7$ for points
1 and 2 respectively upon switching from room to elevated temperatures. The
SHG value integrated over the right hemisphere from the heated dimer was
approximately 2.8 times higher than the signal from the room temperature
dimer. While at low-temperature the SH pattern was almost symmetric because of
the main contribution to the dimer SH resonance from the included
"on"-cylinder, it becomes strongly asymmetric at a higher temperature  due to
the influence of the commensurate radiations from both dimer elements.

Finally, we theoretically proposed an approach to control the directivity of
SH emission from nanoantenna structure with help of thermooptical effect,
which efficiency can be enhanced under proper designing of subwavelength
high-Q systems.

Figure 5:  Second harmonic generation from individual dimer components. The
dependence of SHG intensity normalized to the square of the incident intensity
$\mathrm{I_{SHG}/I_{0}^{2}}$, where
$\mathrm{I_{0}}=10^{13}\>\mathrm{W/m^{2}}$, at two temperature values
separately from the cylinder (solid lines) and prism (dotted lines) that make
up the dimer. The blue lines corresponds to 25o C and the orange line is for
the case of the elevated temperature of $325^{o}$ C. The fundamental
wavelength is 1565 nm and is marked with the vertical black line.

Acknowledgments The work was supported by the Grant Of President of Russian
Federation (grant no. MK-2360.2020.2) and Russian Foundation of Basic Research
(grant no. 20-32-90238). M.P. acknowledges support from the federal academic
leadership program Priority 2030.

Disclosures The authors declare no conflicts of interest.

## References

  * [1] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Phys. Rev. Lett. 7, 118 (1961).
  * [2] M. Yamada, N. Nada, M. Saitoh, and K. Watanabe, Appl. Phys. Lett. 62, 435 (1993).
  * [3] P. Campagnola, Anal. Chem. 83, 3224 (2011).
  * [4] P. J. Campagnola and L. M. Loew, Nat. Biotechnol. 21, 1356 (2003).
  * [5] M. M. Fejer, G. A. Magel, D. H. Jundt, and R. L. Byer, IEEE J. Quantum Electron. 28, 2631 (1992).
  * [6] M. Celebrano, X. Wu, M. Baselli, S. Großmann, P. Biagioni, A. Locatelli, C. De Angelis, G. Cerullo, R. Osellame, B. Hecht, L. Duò, F. Ciccacci, and M. Finazzi, Nat. Nanotechnol. 10, 412 (2015).
  * [7] J. Butet, P.-F. Brevet, and O. J. F. Martin, ACS Nano 9, 10545 (2015).
  * [8] K. Mochizuki, M. Sugiura, H. Yogo, S. Lundgaard, J. Hu, S. H. Ng, Y. Nishijima, S. Juodkazis, and A. Sugita, Micromachines 11, 848 (2020).
  * [9] M. Finazzi, P. Biagioni, M. Celebrano, and L. Duò, Phys. Rev. B 76, 125414 (2007).
  * [10] A. Noor, A. R. Damodaran, I.-H. Lee, S. A. Maier, S.-H. Oh, and C. Ciracì, ACS Photonics 7, 3333 (2020).
  * [11] S. Elkabetz, K. N. Reddy, P. Y. Chen, A. I. Fernández-Domínguez, and Y. Sivan, Phys. Rev. B 103, 075411 (2021).
  * [12] K. Nadolski, E. Benichou, N. Tarnowicz-Staniak, A. Żak, C. Jonin, K. Matczyszyn, and P.-F. Brevet, J. Phys. Chem. C 124, 14797 (2020).
  * [13] E. Gürdal, A. Horneber, A. J. Meixner, D. P. Kern, D. Zhang, and M. Fleischer, Nanoscale 12, 23105 (2020).
  * [14] S. V. Makarov, M. I. Petrov, U. Zywietz, V. Milichko, D. Zuev, N. Lopanitsyna, A. Kuksin, I. Mukhin, G. Zograf, E. Ubyivovk, D. A. Smirnova, S. Starikov, B. N. Chichkov, and Y. S. Kivshar, Nano Lett. 17, 3047 (2017).
  * [15] D. Smirnova and Y. S. Kivshar, Optica 3, 1241 (2016).
  * [16] H. Jiang, Y. Cai, and Z. Han, J. Opt. Soc. Am. B, JOSAB 37, 3146 (2020).
  * [17] C. Gigli, T. Wu, G. Marino, A. Borne, G. Leo, and P. Lalanne, ACS Photonics 7, 1197 (2020).
  * [18] L. Xu, M. Rahmani, D. Smirnova, K. Zangeneh Kamali, G. Zhang, D. Neshev, and A. E. Miroshnichenko, Photonics 5, 29 (2018).
  * [19] T. Zeng, L. Guo, L. Xu, and T. Han, J. Opt. Soc. Am. B, JOSAB 37, 868 (2020).
  * [20] I. Volkovskaya, L. Xu, L. Huang, A. I. Smirnov, A. E. Miroshnichenko, and D. Smirnova, Nanophotonics 9, 3953 (2020).
  * [21] J. Cambiasso, G. Grinblat, Y. Li, A. Rakovich, E. Cortés, and S. A. Maier, Nano Lett. 17, 1219 (2017).
  * [22] M. K. Kroychuk, A. S. Shorokhov, D. F. Yagudin, D. A. Shilkin, D. A. Smirnova, I. Volkovskaya, M. R. Shcherbakov, G. Shvets, and A. A. Fedyanin, Nano Lett. 20, 3471 (2020).
  * [23] G. Grinblat, M. Rahmani, E. Cortés, M. Caldarola, D. Comedi, S. A. Maier, and A. V. Bragas, Nano Lett. 14, 6660 (2014).
  * [24] F. Timpu, N. R. Hendricks, M. Petrov, S. Ni, C. Renaut, H. Wolf, L. Isa, Y. Kivshar, and R. Grange, Nano Lett. 17, 5381 (2017).
  * [25] Z.-J. Yang, Q. Zhao, Y.-H. Deng, D. Zhang, and J. He, Opt. Express 26, 5835 (2018).
  * [26] S. Chervinskii, K. Koskinen, S. Scherbak, M. Kauranen, and A. Lipovskii, Phys. Rev. Lett. 120, 113902 (2018).
  * [27] S. Scherbak and A. A. Lipovskii, J. Phys. Chem. C 122, 15635 (2018).
  * [28] H.-G. Hu, M. Yang, P. Yue, Y.-T. Bai, W.-J. Wang, and S.-D. Liu, J. Phys. D: Appl. Phys. 53, 215101 (2020).
  * [29] H. Linnenbank, Y. Grynko, J. Förstner, and S. Linden, Light Sci. Appl. 5, e16013 (2016).
  * [30] G. Wang, X. Marie, I. Gerber, T. Amand, D. Lagarde, L. Bouet, M. Vidal, A. Balocchi, and B. Urbaszek, Physical review letters 114, 097403 (2015).
  * [31] S. Busschaert, R. Reimann, M. Cavigelli, R. Khelifa, A. Jain, and L. Novotny, ACS Photonics 7, 2482 (2020).
  * [32] R. W. Boyd, _Nonlinear optics_ (Academic press, 2020).
  * [33] L. Kang, H. Bao, and D. H. Werner, Optics Letters 46, 633 (2021).
  * [34] T. Ning, X. Li, Z. Zhang, Y. Huo, Q. Yue, L. Zhao, and Y. Gao, Optics Express 29, 17286 (2021).
  * [35] Z. Han, F. Ding, Y. Cai, and U. Levy, Nanophotonics 10, 1189 (2021).
  * [36] K. Koshelev, S. Kruk, E. Melik-Gaykazyan, J.-H. Choi, A. Bogdanov, H.-G. Park, and Y. Kivshar, Science 367, 288 (2020).
  * [37] G. P. Zograf, M. I. Petrov, D. A. Zuev, P. A. Dmitriev, V. A. Milichko, S. V. Makarov, and P. A. Belov, Nano Lett. 17, 2945 (2017).
  * [38] G. P. Zograf, M. I. Petrov, S. V. Makarov, and Y. S. Kivshar, Advances in Optics and Photonics 13, 643 (2021).
  * [39] P. Yazdanfar, H. Heydarian, and B. Rashidian, J. Opt. Soc. Am. B, JOSAB 37, 2238 (2020).
  * [40] T. V. Tsoulos and G. Tagliabue, Nanophotonics 9, 3849 (2020).
  * [41] Y.-S. Duh, Y. Nagasaki, Y.-L. Tang, P.-H. Wu, H.-Y. Cheng, T.-H. Yen, H.-X. Ding, K. Nishida, I. Hotta, J.-H. Yang, Y.-P. Lo, K.-P. Chen, K. Fujita, C.-W. Chang, K.-H. Lin, J. Takahara, and S.-W. Chu, Nat. Commun. 11, 1 (2020).
  * [42] M. Rahmani, L. Xu, A. E. Miroshnichenko, A. Komar, R. Camacho-Morales, H. Chen, Y. Zárate, S. Kruk, G. Zhang, D. N. Neshev, and Y. S. Kivshar, Adv. Funct. Mater. 27, 1700580 (2017).
  * [43] K. Z. Kamali, L. Xu, J. Ward, K. Wang, G. Li, A. E. Miroshnichenko, D. Neshev, and M. Rahmani, Small 15, 1805142 (2019).
  * [44] T. Zhang, Y. Che, K. Chen, J. Xu, Y. Xu, T. Wen, G. Lu, X. Liu, B. Wang, X. Xu, Y.-S. Duh, Y.-L. Tang, J. Han, Y. Cao, B.-O. Guan, S.-W. Chu, and X. Li, Nat. Commun. 11 (2020).
  * [45] M. Celebrano, D. Rocco, M. Gandolfi, A. Zilli, F. Rusconi, A. Tognazzi, A. Mazzanti, L. Ghirardini, E. A. A. Pogna, L. Carletti, C. Baratto, G. Marino, C. Gigli, P. Biagioni, L. Duò, G. Cerullo, G. Leo, G. Della Valle, M. Finazzi, and C. De Angelis, Optics Letters 46, 2453 (2021).
  * [46] L. Carletti, A. Locatelli, D. Neshev, and C. De Angelis, ACS Photonics 3, 1500 (2016).
  * [47] C. Gigli, G. Marino, A. Artioli, D. Rocco, D. Rocco, C. De Angelis, C. De Angelis, J. Claudon, J.-M. Gérard, and G. Leo, Optica 8, 269 (2021).
  * [48] I. Shoji, T. Kondo, A. Kitamoto, M. Shirane, and R. Ito, J. Opt. Soc. Am. B, JOSAB 14, 2268 (1997).
  * [49] M. Ohashi, T. Kondo, R. Ito, S. Fukatsu, Y. Shiraki, K. Kumata, and S. S. Kano, J. Appl. Phys. 74, 596 (1993).
  * [50] D. Aspnes, S. Kelso, R. Logan, and R. Bhat, Journal of applied physics 60, 754 (1986).
  * [51] S. Gehrsitz, F. K. Reinhart, C. Gourgon, N. Herres, A. Vonlanthen, and H. Sigg, J. Appl. Phys. 87, 7825 (2000).
  * [52] S. Adachi, J. Appl. Phys. 58, R1 (1985).
  * [53] M. A. Afromowitz, Journal of Applied Physics 44, 1292 (1973).
  * [54] K. Frizyuk, I. Volkovskaya, D. Smirnova, A. Poddubny, and M. Petrov, Physical Review B 99, 075425 (2019).
  * [55] Y. Zeng, H. Qian, M. J. Rozin, Z. Liu, and A. R. Tao, Adv. Funct. Mater. 28, 1803019 (2018).
  * [56] K. Thyagarajan, S. Rivier, A. Lovera, and O. J. F. Martin, Opt. Express 20, 12860 (2012).
  * [57] S. Park, J. W. Hahn, and J. Y. Lee, Opt. Express 20, 4856 (2012).
  * [58] J. Wang, M. Clementi, M. Minkov, A. Barone, J.-F. Carlin, N. Grandjean, D. Gerace, S. Fan, M. Galli, and R. Houdré, Optica 7, 1126 (2020).
  * [59] S. Gladyshev, K. Frizyuk, and A. Bogdanov, Physical Review B 102, 75103 (2020).
  * [60] S. Liu, M. B. Sinclair, S. Saravi, G. A. Keeler, Y. Yang, J. Reno, G. M. Peake, F. Setzpfandt, I. Staude, T. Pertsch _et al._ , Nano letters 16, 5426 (2016).
  * [61] M. Aouassa, E. Mitsai, S. Syubaev, D. Pavlov, A. Zhizhchenko, I. Jadli, L. Hassayoun, G. Zograf, S. Makarov, and A. Kuchmizhak, Applied Physics Letters 111, 243103 (2017).
  * [62] A. E. Miroshnichenko and M. I. Tribelsky, Physical review letters 120, 033902 (2018).
  * [63] O. Pashina, D. Ryabov, G. Zograf, S. Makarov, and M. Petrov, “Ultrafast laser heating of non-plasmonic nanocylinders,” in _Journal of Physics: Conference Series,_ , vol. 2015 (IOP Publishing, 2021), p. 012104.
  * [64] H. Yao, P. G. Snyder, and J. A. Woollam, J. Appl. Phys. 70, 3261 (1991).
  * [65] H. Yao, P. G. Snyder, K. Stair, and T. Bird, MRS Online Proceedings Library (OPL) 242 (1992).
  * [66] G. Baffou, _Thermoplasmonics: Heating metal nanoparticles using light_ (Cambridge University Press, 2017).
  * [67] K. Frizyuk, JOSA B 36, F32 (2019).