# Demonstration of a plasmonic nonlinear pseudo-diode

Sergejs Boroviks<EMAIL_ADDRESS>Nanophotonics and Metrology
Laboratory, Swiss Federal Institute of Technology Lausanne (EPFL), Lausanne,
Switzerland  Andrei Kiselev  Nanophotonics and Metrology Laboratory, Swiss
Federal Institute of Technology Lausanne (EPFL), Lausanne, Switzerland  Karim
Achouri  Nanophotonics and Metrology Laboratory, Swiss Federal Institute of
Technology Lausanne (EPFL), Lausanne, Switzerland  Olivier J.F. Martin
<EMAIL_ADDRESS>Nanophotonics and Metrology Laboratory, Swiss Federal
Institute of Technology Lausanne (EPFL), Lausanne, Switzerland

(January 30, 2023)

###### Abstract

We demonstrate a nonlinear plasmonic metasurface that exhibits strongly
asymmetric second-harmonic generation: nonlinear scattering is efficient upon
excitation in one direction and it is substantially suppressed when the
excitation direction is reversed, thus enabling a diode-like functionality. A
significant (approximately $10\text{\,}\mathrm{d}\mathrm{B}$) extinction ratio
of SHG upon opposite excitations is measured experimentally and those findings
are substantiated with full-wave simulations. The combination of two commonly
used metals – aluminium and silver – produces a material composition asymmetry
that results into a bianisotropic response of the system, as confirmed by
performing homogenization analysis and extracting an effective susceptibility
tensor. Finally, we discuss the implications of our results from the more
fundamental perspectives of reciprocity and time-reversal asymmetry.

Plasmonics, Metasrufaces, Bianisotropy, Time-reversal asymmetry, Nonlocality,
Second-Harmonic Generation

High-performance nanoscale devices that allow transmission of light only in
one direction – optical isolators – remain a long-coveted research objective
for optical engineers. This problem is nontrivial due to the fundamental
property of electromagnetic waves: in linear time-invariant (LTI) media and in
the absence of an external time-odd bias, such as a magnetic field, they
propagate reciprocally, i.e. the same way in the forward and backward
directions. This property is linked with the time-reversal symmetry of the
macroscopic Maxwell’s equations and can be shown via the Lorentz reciprocity
theorem, which specifically applies to LTI media Caloz  _et al._ (2018);
Asadchy  _et al._ (2020); Achouri and Martin (2021). However, despite recent
comprehensive publications on this topic Jalas  _et al._ (2013); Sounas and
Alu (2017); Caloz  _et al._ (2018); Asadchy  _et al._ (2020); Sigwarth and
Miniatura (2022), there remains a tangible confusion in the community about
the difference between true nonreciprocity and deceptively similar time-
reversal asymmetric response. For example, time-invariant and bias-less lossy
systems may exhibit contrast upon excitation from opposite directions, but
they do not qualify as optical isolators since they possess a symmetric
scattering matrix and thus obey Lorentz reciprocity Fan  _et al._ (2012).
Furthermore, in the case of devices based on nonlinear effects, the
distinction between true and pseudo-isolators is even more intricate. In
particular, devices based on Kerr-type nonlinearities Cotrufo  _et al._ (2021)
are intrinsically limited by dynamic reciprocity: they can only perform as
pseudo-isolators, since they do not exhibit unidirectional transmission upon
simultaneous excitation from opposite directions Shi  _et al._ (2015);
Fernandes and Silveirinha (2018). One aim of this work is to explore
possibilities to overcome this limitation and demonstrate how it can be turned
into an advantage with an appropriate application.

In that context, photonic metasurfaces – artificial planar materials
constituted of subwavelength elements – have been identified as a promising
platform for the realization of miniature optical isolators or asymmetric
devices Shaltout  _et al._ (2019). To this end, let us highlight recent
progress in the development of two classes of metasurfaces – nonlinear and
bianisotropic metasurfaces. These two classes are particularly relevant to the
scope of our work, since combining their features enables realization of
unconventional functionalities, such as aforementioned nonlinearly induced
nonreciprocity Menzel  _et al._ (2010); Mahmoud  _et al._ (2015); Lawrence
_et al._ (2018); Chen  _et al._ (2021); Cheng  _et al._ (2021), directional
harmonic generation Yang  _et al._ (2017a); Xu  _et al._ (2020); Nauman  _et
al._ (2021) and nonlinear beam shaping Tymchenko  _et al._ (2016); Bar-David
and Levy (2019).

Nonlinear metasurfaces Minovich  _et al._ (2015); Li  _et al._ (2017); Krasnok
_et al._ (2018) have the potential to replace bulky optical crystals and thus
miniaturize nonlinear optical devices. Among other applications, plasmonic
metasurfaces have proven to be interesting for second-harmonic generation
(SHG) Lee  _et al._ (2014); Kauranen and Zayats (2012); Butet  _et al._
(2015), which is a second-order nonlinear optical process in which an
excitation wave with frequency $\omega$ is converted into a wave with double
frequency $2\omega$ Boyd (2020). However, the second-order nonlinear response
of plasmonic metals is weak due to their centrosymmetric crystal structure,
which is only broken at the surface, giving rise to a non-vanishing surface
normal component of the second-order susceptibility tensor
$\chi_{\perp\perp\perp}^{(2)}$. Yet, the overall SHG efficiency remains small
due to the reduced interaction volume: essentially, the nonlinear process
occurs within the few atomic layers at the metal surface, since the bulk metal
is opaque for visible and infrared light and its bulk second-order response is
vanishing. Nevertheless, this limitation can be partially overcome by the
virtue by virtue of the field enhancement associated with surface plasmon
resonances at metal surfaces. Thus, various SHG enhancement schemes were
proposed for plasmonic metasurfaces, based on multipolar resonances
Chandrasekar  _et al._ (2015); Kruk  _et al._ (2015); Smirnova and Kivshar
(2016); Bernasconi  _et al._ (2016); Butet  _et al._ (2017); Yang  _et al._
(2017b); Kiselev  _et al._ (2019), plasmonic lattice resonances Gupta  _et
al._ (2021); Abir  _et al._ (2022) and even light-induced centrosymmetry
breaking Li  _et al._ (2021).

On the other hand, bianisotropic metasurfaces allow engineering the
polarization response to realize highly efficient refraction devices through
the combination of electric and magnetic effects Asadchy  _et al._ (2018);
Achouri and Caloz (2021). The bianisotropic response, which emerges in
structures with broken spatial symmetries Achouri  _et al._ (2022a), implies
that the material acquires magnetic polarization upon excitation with an
electric field, and vice versa, electric polarization is produced by a
magnetic field. Such a magneto-electric coupling gives rise to the spatial
dispersion (i.e. wavevector-dependent response) that enables an excitation
angle-dependent operation Overvig and Alù (2022). For example, in lossy
systems, it may lead to asymmetric reflection and absorption, which will be
discussed further in relation to our work.

Figure 1: Comparison of conventional and asymmetric SHG: (a) symmetric SHG
from a nonlinear (NL) crystal; (b) asymmetric SHG from a nonlinear
bianisotropic metasurface.

In this work, we demonstrate theoretically and experimentally a plasmonic
metasurface that exhibits asymmetric SHG. The operation of the device is
conceptually depicted in Fig. 1: in contrast to a conventional nonlinear
crystal, second-harmonic (SH) is efficiently generated only upon one
excitation direction, which essentially, enables a nonlinear optical pseudo-
diode functionality (to be distinguished from optical isolators and pseudo-
isolators). Such an asymmetric response imposes a structural asymmetry of the
system and previously proposed theoretical designs with similar
functionalities have relied on a geometric asymmetry, which might be difficult
to realize experimentally Poutrina and Urbas (2016); Mobini  _et al._ (2021);
Jin and Argyropoulos (2020); Kim (2021); Liu  _et al._ (2021). Here, we take a
different route and implement a structural asymmetry through the utilization
of two common plasmonic meterials – silver (Ag) and aluminium (Al) – in a
metasurface and show that substantial direction-dependent SHG (up to approx.
$16.9\text{\,}\mathrm{d}\mathrm{B}$ in theory and approx.
$10\text{\,}\mathrm{d}\mathrm{B}$ in experiment). A major advantage of this
two-dimensional design is that such a material asymmetry is relatively easy to
implement using standard nanofabrication techniques, e.g. single-exposure
electron-beam lithography (EBL) Abasahl  _et al._ (2021). Furthermore, the
combination of plasmonic metals is known to enhance nonlinear processes Wang
_et al._ (2019, 2021). To the best of our knowledge, this is the first
experimental demonstration of a plasmonic metasurface for asymmetric SHG,
although we note that in a recent experimental demonstration Kruk et al.
utilized a combination of dielectric nonlinear materials for third-harmonic
generation Kruk  _et al._ (2022). Additionally, we perform homogenization
analysis of the metasurface to extract effective susceptibilities and reveal
bianiostropic property of our metasurface. Finally, we discuss the fundamental
implications of our results in the context of nonreciprocity.

The building block of the metasurface – the meta-atom – is schematically
depicted in Fig. 2a. It is comprised of two T-shaped nanostructures made of Al
and Ag that are stacked one on top of the other and separated by a thin
silicone dioxide (SiO2) spacer. These nanostructures are embedded in SiO2 and
arranged in a square lattice with the period of
$\Lambda=$$250\text{\,}\mathrm{nm}$. Such a periodicity is sufficiently small
to avoid diffraction in both linear and nonlinear regimes, as the metasurface
is designed for the excitation with the vacuum wavelength of
$\lambda_{0}=$$800\text{\,}\mathrm{nm}$ (the effective wavelength in SiO2 is
$\sim$$537\text{\,}\mathrm{nm}$) and SHG at
$\lambda_{\textrm{SH}}=$$400\text{\,}\mathrm{nm}$
($\sim$$268\text{\,}\mathrm{nm}$ in SiO2).

As shown in Fig. 2b, we consider two different excitation conditions that are
indicated with red thick arrows: forward (in the direction along the
$+z$-axis) and backward (along the $-z$-axis) propagating plane waves that are
$x$-polarized. In the linear regime, each of the two waves gives rise to
transmitted (red solid arrows) and reflected (red dashed arrows) waves, which
are labeled as forward-excited reflection (FR) and transmission (FT), or
backward-excited reflection (BR) and transmission (BT). Additionally, both
excitations produce signals at the SH frequency (shown with blue arrows). For
the SH signals, we use the same naming convention as the waves produced by
linear scattering, Fig. 2b. For the reflected and transmitted waves at the
excitation frequency, we measure the co-polarized $x$-component of the
electric field, whereas for the SHG waves, the cross-polarized $y$-component
is measured, as it is found to be dominant (see Fig. S3 in the Supporting
Information).

T-shaped meta-atoms provide almost independent control of the spectral
positions for the resonances both at the excitation and SH frequencies by
varying the lateral dimensions $L_{x}$ and $L_{y}$ Czaplicki  _et al._ (2015).
As can be seen from Fig. S1 in the Supporting Information, for a fixed
wavelength, the transmission in the linear regime is tuned by varying $L_{x}$.
In the nonlinear regime, the transmission and reflection are controlled by
both $L_{x}$ and $L_{y}$. Importantly, for forward excitation, the maximum in
SHG transmission coincides with the minimum in linear transmission (compare
panels a and b in Fig. S1 in the Supporting Information). The other geometric
parameters $L_{\textrm{s}}$, $D$, $t_{\ce{Ag}}$ and $t_{\ce{Al}}$ do not have
a strong influence on the resonance wavelength of the fundamental mode,
however they affect the scattering cross-section of the meta-atoms via the
retardation effects Kottmann and Martin (2001), which, in turn determines the
overall transmission and SHG intensity (see Fig. S2 in the Supporting
Information). The sidewalls of the meta-atom are tilted by
$10\text{\,}\mathrm{\SIUnitSymbolDegree}$ and the edges and corners are
rounded with a $5\text{\,}\mathrm{nm}$ radius to mimic the experimentally
fabricated structures, as discussed below.

We select $L_{x}=$$135\text{\,}\mathrm{nm}$,
$L_{y}=$$195\text{\,}\mathrm{nm}$, $L_{\textrm{s}}=$$25\text{\,}\mathrm{nm}$
and $D=t_{\ce{Ag}}=t_{\ce{Al}}=$$50\text{\,}\mathrm{nm}$, since these
parameters maximize SHG upon forward excitation at the design wavelength. Such
meta-atom dimensions result in minimal transmission in the linear regime and
sufficiently high extinction ratio of SHG upon forward and backward excitation
(see the parametric sweeps in Fig. S1 in the Supporting Information).
Furthermore, in the $L_{x}$ and $L_{y}$ parameter space, the forward-
excitation SHG peak is broad, which implies that the metasurface efficiency is
weakly sensitive to deviations from the nominal dimensions, thus easing
nanofabrication tolerances.

The simulations are performed in two steps using a custom-developed numerical
electromagnetic solver based on the surface integral equation Gallinet  _et
al._ (2010); Butet  _et al._ (2013). First, the linear fields are computed
with a plane-wave excitation and periodic boundary conditions. For the SHG
simulations, the nonlinear surface polarization
$P_{\perp}^{(2\omega)}=\chi^{(2)}_{\perp\perp\perp}E_{\perp}^{\omega}E_{\perp}^{\omega}$
is used as a source, where the normal components of the surface fields
$E_{\perp}^{\omega}$ are obtained from the linear simulations.

The simulated reflectance and transmittance in the linear and SHG regimes are
shown in Fig. 2c and d. In the simulations, we use interpolated values
$\varepsilon_{\ce{Al}}$ and $\varepsilon_{\ce{Ag}}$ of the experimental
permittivity data from McPeak et al. McPeak  _et al._ (2015), and for the
permittivity of the background medium we use $\varepsilon_{\ce{SiO2}}=2.22$.
Among the noble metals, Ag is known to have the lowest losses at optical
frequencies, whereas Al has recently attracted attention as a low cost
alternative plasmonic material Castro-Lopez  _et al._ (2011); Knight  _et al._
(2014); Gérard and Gray (2014); Thyagarajan  _et al._ (2016). Apart from its
low cost, Al is known to have the highest second-order nonlinear
susceptibility among the plasmonic materials, in particular its surface normal
component $\chi^{(2)}_{\perp\perp\perp}$ Krause  _et al._ (2004), it also
exhibits an interband transition-related absorption peak at
$800\text{\,}\mathrm{nm}$ (see Fig. S4 in the Supporting Information).

Figure 2: Design and simulated performance of the nonlinear bianisotropic
metasurface. (a) Schematics of the system and the considered forward-
excitation transmission (FT) and reflection (FR), as well as backward-
excitation transmission (BT) and reflection (BR); thick solid red arrows
indicate the excitation waves; thin solid (dashed) arrows indicate the
transmitted (reflected) waves at the excitation frequency in red and at the SH
frequency in blue. (b) Schematic drawing of the metasurface unit cell in
isometric-, top- and side-views with indicated geometric and material
parameters. Simulated metasurface reflectance and transmittance (c) in the
linear regime and (d) at the SH frequency. Relevant components of the
extracted (e) linear and (f) nonlinear effective susceptibility tensors.

As shown in Fig. 2c, in the linear regime, the transmission $T$ for both
forward and backward excitations is exactly the same, as imposed by
reciprocity. However, the reflection $R$ and absorption $A$, which are related
to transmission as $A+R=1-T$, depend on the excitation direction, as they are
not restricted by reciprocity and depend on the spatial asymmetry of the
system. The asymmetric reflection and absorption of the system can be analyzed
by considering an isolated meta-atom. As can be seen in Fig. S5c and d in the
Supporting Information, forward and backward excitations give rise to two
distinct electric field distributions. In particular, the electric field
concentration in the Al part of the structure is strongly dependent on the
excitation direction. Although the response is primarily dipolar for both
excitations (see Fig. S6a and b in the Supporting Information), this results
in asymmetric linear scattering and absorption cross-sections, which is a
characteristic of bianisotropic systems Cheng  _et al._ (2021). In fact, it is
presence of the losses that enables asymmetric scattering when the structure
is illuminated from opposite directions, whereas the extinction cross-section
remains exactly the same, as imposed by reciprocity Sounas and Alù (2014).

In turn, the SHG response that is plotted in Fig. 2d, has an even stronger
dependence on the excitation direction: both nonlinear FT and RT are more than
two orders of magnitude stronger than the BT an BR at
$400\text{\,}\mathrm{nm}$. A multipolar analysis of an isolated meta-atom (see
Fig. S6c and d in the Supporting Information), shows that the electric dipolar
and quadruplar modes are excited more efficiently at $400\text{\,}\mathrm{nm}$
upon forward excitation. This is due to the aforementioned different electric-
field distributions at the surface of the T-shaped particles, that become the
sources for the SHG.

To further elucidate the significance of bianisotropy in such an asymmetric
response, we extracted the effective susceptibilities from the simulated
electromagnetic fields following the previously documented procedure of
metasurface homogenization analysis Achouri  _et al._ (2017, 2018, 2022b).
Briefly, the expressions for nonlinear susceptibilities are derived from the
generalized sheet transition conditions and are calculated using the simulated
reflected and transmitted fields upon different excitation conditions at
$\omega$ and $2\omega$ frequencies.

In Fig. 2e and f we plot the extracted effective susceptibility tensor
elements that are relevant to the considered excitation conditions. For both
linear and nonlinear susceptibilities, the magneto-electric coupling
(corresponding to the terms with mixed “e” and “m” subscripts in Fig. 2e and
f) is non-negligible. The asymmetric response becomes apparent by noting that
the induced linear and nonlinear polarizations are given by

$\displaystyle\mathbf{P}^{\omega}=\overline{\overline{\mathbf{\chi}}}_{\textrm{ee}}^{\>\omega}\cdot\mathbf{E}^{\omega}+\overline{\overline{\mathbf{\chi}}}_{\textrm{em}}^{\>\omega}\cdot\mathbf{H}^{\omega},$
(1a)
$\displaystyle\mathbf{P}^{2\omega}=\overline{\overline{\mathbf{\chi}}}_{\textrm{ee}}^{\>2\omega}\cdot\mathbf{E}^{2\omega}+\overline{\overline{\mathbf{\chi}}}_{\textrm{em}}^{\>2\omega}\cdot\mathbf{H}^{2\omega}+\overline{\overline{\mathbf{\chi}}}_{\textrm{eee}}^{\>\omega}:\mathbf{E}^{\omega}\mathbf{E}^{\omega}+\overline{\overline{\mathbf{\chi}}}_{\textrm{eem}}^{\>\omega}:\mathbf{E}^{\omega}\mathbf{H}^{\omega}+\overline{\overline{\mathbf{\chi}}}_{\textrm{emm}}^{\>\omega}:\mathbf{H}^{\omega}\mathbf{H}^{\omega}.$
(1b)

In the linear regime, the non-negligible magneto-electric coupling term
$\chi_{\textrm{me}}$ results in an asymmetric absorption and reflection. As
for the nonlinear effective susceptibility tensors, the dominant components
are $\chi_{\textrm{mem}}^{yyx}$ and $\chi_{\text{eem}}^{xyx}$, which relate
magnetic/electric excitations with electric/magnetic responses along
orthogonal directions and result in strongly asymmetric SHG.

To verify experimentally this asymmetric nonlinear response, we fabricated and
characterized a metasurface device. Instead of the widespread lift-off
process, we employ the ion beam etching (IBE) technique which enables the
fabrication of stratified nanostructures, in particular metal-dielectric
composites, with sharper features Ray  _et al._ (2020); Abasahl  _et al._
(2021). The schematic flowchart of the fabrication process is shown in Fig.
3a. We use a $150\text{\,}\mathrm{\SIUnitSymbolMicro m}$-thick D 263 glass
wafer (Schott) which is coated with $50\text{\,}\mathrm{nm}$-thick Al and
$25\text{\,}\mathrm{nm}$ SiO2 films using RF sputtering (Pfeiffer SPIDER 600).
Next, we deposit a $50\text{\,}\mathrm{nm}$ thick Ag layer using an e-beam
assisted evaporator (Alliance-Concept EVA 760). The T-shaped pattern arrays
are exposed in the hydrogen silsesquioxane (HSQ, XR-1541-006 from DuPont),
which is a negative tone e-beam resist, using electron beam lithography (Raith
EBPG5000+). The formation of the exposed patterns in the thin films is
performed using a low-power argon IBE (Veeco Nexus IBE350, operated at a
$300\text{\,}\mathrm{V}$ acceleration voltage). An important point for this
last step is the pulsed IBE operation: 10 s of etching followed by 30 s of
cooling to avoid damaging the sample by substrate overheating. The typical
overall IBE process time is $160\text{\,}\mathrm{s}$, and the etching depth is
controlled in-situ using a mass-spectrometer, which allows real-time
monitoring of the etched material composition: the etching process is stopped
as soon as the Al flux drops to a minimum. The fabrication results are shown
in the scanning electron microscope (SEM) images in Fig. 3b-d. The morphology
of the fabricated structure can be inspected in Fig. 3c: intrinsically, the
IBE process results in tilted sidewalls (approx.
$10\text{\,}\mathrm{\SIUnitSymbolDegree}$) and rounded corners and edges.
Although such features are typically undesired, they are not expected to
degrade the performance of the metasurface, as these were taken into account
in the simulations. In turn, the layered material composition can be well
identified in the image acquired with the back-scattered electron (BSE)
detector in Fig. 3d. In the last fabrication step, we cover the metallic
nanostructures with a thick SiO2 layer (approx. $300\text{\,}\mathrm{nm}$)
which serves two purposes: it acts as a protective layer preventing
degradation of the Al and Ag nanostructures, and simplifies the physical
conditions by having identical permittivities above and below the metasurface.

Figure 3: Fabrication of the bimetallic metasurface. (a) Flowchart of the
fabrication: 1. Initial substrate; 2. Al, SiO2, Ag and HSQ thin films
deposition; 3. E-beam exposure; 4. IBE; 5. Covering with a thick SiO2 film.
SEM images of the fabricated structure acquired using different detectors and
tilt angles: (b) top view SE (scale bar: $1\text{\,}\mathrm{\SIUnitSymbolMicro
m}$. (c) $45\text{\,}\mathrm{\SIUnitSymbolDegree}$-tilted view SE (scale bar:
$200\text{\,}\mathrm{nm}$. (d)
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$-tilted view BSE (scale bar:
$200\text{\,}\mathrm{nm}$.

The experimental setup and the results for the optical characterization of the
fabricated sample are shown in Fig. 4. As an excitation light source, we use a
mode-locked Ti:Saph laser that outputs approx. $120\text{\,}\mathrm{fs}$
pulses with a central wavelength of $800\text{\,}\mathrm{nm}$. The excitation
light is weakly focused onto the metasurface with a low magnification
objective ($\text{NA}=0.1$), which results in a focal spot with a
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ FWHM mimicking the plane wave
excitation used in the simulations. The spectrum of the nonlinearly generated
light is shown in Fig. 3a. Apart from the characteristic SHG peak at
$400\text{\,}\mathrm{n}\mathrm{m}$, it has a tail at longer wavelengths, which
is attributed to nonlinear photo-luminescence (NPL).

As an interesting side-effect, we note that the NPL signal is substantially
larger for BT than for FT. This fact can be explained by the peculiarity of
the two-photon absorption mechanism in metals that induces the NPL. As opposed
to the coherent nature of two-photon absorption in molecules or dielectrics,
in metals it can be regarded as a cascaded process. Specifically, two photons
are absorbed sequentially rather than simultaneously Beversluis  _et al._
(2003); Mühlschlegel  _et al._ (2005); Biagioni  _et al._ (2009). Absorption
of the first photon gives rise to an intraband transition in the conduction
band and creates a vacancy below the Fermi level. Thus, the second photon
results in an interband transition that fills the vacancy in the conduction
band and creates one in the valence band. Both of these photon absorption
steps are linear, but result in an effective nonlinearity. Thus, higher linear
absorption upon backward excitation (see Fig. S5 and discussion above),
results in a higher probability of two-photon absorption and subsequent NPL,
which is consistent with our observations.

Figure 4: Optical characterization of the metasurface: (a) measurement setup.
(b) Excitation spectrum. (c) Nonlinear spectra.

Such asymmetric behaviour is sometimes referred to as ”nonreciprocal SHG”,
both in the metasurfaces Poutrina and Urbas (2016) and solid state physics
Toyoda  _et al._ (2021); Mund  _et al._ (2021) communities. We share the view
that such a nomenclature is improper in the case of SHG, since the concept of
nonrecipocity is not well-defined for nonlinear optics Trzeciecki and Hübner
(2000); Sounas and Alù (2017); Achouri  _et al._ (2018). For any $N$-port
system, the Lorentz reciprocity implies the symmetry of the scattering matrix
$\overline{\overline{\mathbf{S}}}^{T}=\overline{\overline{\mathbf{S}}}$, where
T denotes the transpose operator. In the case of a two-port system like that
considered in this work in the linear regime, the scattering matrix is given
by

$\overline{\overline{\mathbf{S}}}=\begin{bmatrix}S_{11}&S_{12}\\\
S_{21}&S_{22}\end{bmatrix},$ (2)

and reciprocity requires that the transmission coefficients $S_{12}$ and
$S_{21}$ are equal. However, it does not impose any limitations on the
reflection coefficients $S_{11}$ and $S_{22}$. This is true for our system in
the linear regime, since the transmissions for forward and backward
excitations are equal, while the reflections are asymmetric.

However, in the nonlinear regime, our metasurface cannot be regarded as a two-
port system anymore, since the SH emission represents a distinct
electromagnetic mode. Therefore, this system must be at least considered as a
4-port system (assuming that higher-order harmonic generation is negligible),
represented with the following scattering matrix:

$\overline{\overline{\mathbf{S}}}=\begin{bmatrix}\overline{\overline{\mathbf{S}}}^{\omega{\scriptscriptstyle\rightarrow}\omega}&\overline{\overline{\mathbf{S}}}^{2\omega{\scriptscriptstyle\rightarrow}\omega}\\\
\overline{\overline{\mathbf{S}}}^{\omega{\scriptscriptstyle\rightarrow}2\omega}&\overline{\overline{\mathbf{S}}}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}\end{bmatrix}\\\
=\begin{bmatrix}S_{11}^{\omega{\scriptscriptstyle\rightarrow}\omega}&S_{12}^{\omega{\scriptscriptstyle\rightarrow}\omega}&S_{11}^{2\omega{\scriptscriptstyle\rightarrow}\omega}&S_{12}^{2\omega{\scriptscriptstyle\rightarrow}\omega}\\\
S_{21}^{\omega{\scriptscriptstyle\rightarrow}\omega}&S_{22}^{\omega{\scriptscriptstyle\rightarrow}\omega}&S_{21}^{2\omega{\scriptscriptstyle\rightarrow}\omega}&S_{22}^{2\omega{\scriptscriptstyle\rightarrow}\omega}\\\
S_{11}^{\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{12}^{\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{11}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{12}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}\\\
S_{21}^{\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{22}^{\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{21}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}&S_{22}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}\end{bmatrix},$
(3)

which describes both linear transmission/reflection at frequencies $\omega$
and $2\omega$, as well as nonlinear processes
$\omega{\scriptscriptstyle\rightarrow}2\omega$ and
$2\omega{\scriptscriptstyle\rightarrow}\omega$.

In our experiment, we do not directly probe
$S_{21}^{\omega{\scriptscriptstyle\rightarrow}2\omega}\stackrel{{\scriptstyle?}}{{=}}S_{12}^{2\omega{\scriptscriptstyle\rightarrow}\omega}$,
where $S_{12}^{2\omega{\scriptscriptstyle\rightarrow}\omega}$ parameters
corresponds to the excitation at SH frequency and generation of a wave at
frequency $\omega$. In fact, this process is known as known as parametric down
conversion and it has an extremely low efficiency in comparison with SHG Boyd
(2020). Probing this equality, as well as equality of 8 other parameters that
are flipped by the $\overline{\overline{\mathbf{S}}}^{T}$ operation, namely
$S_{21}^{\omega{\scriptscriptstyle\rightarrow}\omega}\stackrel{{\scriptstyle?}}{{=}}S_{12}^{\omega{\scriptscriptstyle\rightarrow}\omega}$,
$S_{11}^{\omega{\scriptscriptstyle\rightarrow}2\omega}\stackrel{{\scriptstyle?}}{{=}}S_{11}^{2\omega{\scriptscriptstyle\rightarrow}\omega}$,
$S_{12}^{\omega{\scriptscriptstyle\rightarrow}2\omega}\stackrel{{\scriptstyle?}}{{=}}S_{21}^{2\omega{\scriptscriptstyle\rightarrow}\omega}$,
$S_{22}^{\omega{\scriptscriptstyle\rightarrow}2\omega}\stackrel{{\scriptstyle?}}{{=}}S_{22}^{2\omega{\scriptscriptstyle\rightarrow}\omega}$
and
$S_{21}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}\stackrel{{\scriptstyle?}}{{=}}S_{12}^{2\omega{\scriptscriptstyle\rightarrow}2\omega}$
stand for a true reciprocity test in a four-port system. Instead, within our
experiment we show that
$S_{21}^{\omega{\scriptscriptstyle\rightarrow}2\omega}\neq
S_{12}^{\omega{\scriptscriptstyle\rightarrow}2\omega}$, which corresponds to
an asymmetric nonlinear scattering process that is reciprocal. Yet, a rigorous
probing of reciprocity in a nonlinear system would require sophisticated
experiments that involve simultaneous excitation with the two waves at
frequencies $\omega$ and $2\omega$ and precise control over their amplitude
and phase Trzeciecki and Hübner (2000). Nevertheless, we assert that our
device essentially functions as a nonlinear optical pseudo-diode, allowing the
transmission of SH signal only in one direction, which is a desired
functionality for various signal processing applications Willner  _et al._
(2014).

In summary, we have demonstrated that strongly asymmetric SHG can be achieved
in a plasmonic metasurface that is comprised of two common plasmonic metals –
aluminium and silver. The structural asymmetry created by the material
contrast results in a strong dependence on the excitation direction, with an
extinction ratio of approx. $16.9\text{\,}\mathrm{d}\mathrm{B}$ in theory and
approx. $10\text{\,}\mathrm{d}\mathrm{B}$ in the experiment. We anticipate
that our findings can pave the way for further developments in the field of
nonlinear bianisotropic and nonreciprocal devices, as well as inspire novel
plasmonic devices with unrivaled functionalities.

## acknowledgement

The authors thank Christian Santschi and Zdenek Benes for their valuable
advises on nanofabrication.
Funding from the Swiss National Science Foundation (grant PZ00P2_193221) is
gratefully acknowledged.

## References

  * Caloz  _et al._ (2018) C. Caloz, A. Alù, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Léger, Phys. Rev. Applied 10, 047001 (2018).
  * Asadchy  _et al._ (2020) V. S. Asadchy, M. S. Mirmoosa, A. Díaz-Rubio, S. Fan, and S. A. Tretyakov, Proceedings of the IEEE 108, 1684 (2020).
  * Achouri and Martin (2021) K. Achouri and O. J. F. Martin, Phys. Rev. B 104, 165426 (2021).
  * Jalas  _et al._ (2013) D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popovic, A. Melloni, J. D. Joannopoulos, M. Vanwolleghem, C. R. Doerr, and H. Renner, Nature Photonics 7, 579 (2013).
  * Sounas and Alu (2017) D. L. Sounas and A. Alu, Nature Photonics 11, 774 (2017).
  * Sigwarth and Miniatura (2022) O. Sigwarth and C. Miniatura, AAPPS Bulletin 32 (2022).
  * Fan  _et al._ (2012) S. Fan, R. Baets, A. Petrov, Z. Yu, J. D. Joannopoulos, W. Freude, A. Melloni, M. Popović, M. Vanwolleghem, D. Jalas, M. Eich, M. Krause, H. Renner, E. Brinkmeyer, and C. R. Doerr, Science 335, 38 (2012).
  * Cotrufo  _et al._ (2021) M. Cotrufo, S. A. Mann, H. Moussa, and A. Alù, IEEE Transactions on Microwave Theory and Techniques 69, 3569 (2021).
  * Shi  _et al._ (2015) Y. Shi, Z. Yu, and S. Fan, Nature Photonics 9, 388 (2015).
  * Fernandes and Silveirinha (2018) D. E. Fernandes and M. G. Silveirinha, IEEE Antennas and Wireless Propagation Letters 17, 1953 (2018).
  * Shaltout  _et al._ (2019) A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, Science 364, eaat3100 (2019).
  * Menzel  _et al._ (2010) C. Menzel, C. Helgert, C. Rockstuhl, E.-B. Kley, A. Tünnermann, T. Pertsch, and F. Lederer, Phys. Rev. Lett. 104, 253902 (2010).
  * Mahmoud  _et al._ (2015) A. M. Mahmoud, A. R. Davoyan, and N. Engheta, Nature Communications 6, 1 (2015).
  * Lawrence  _et al._ (2018) M. Lawrence, D. R. Barton, and J. A. Dionne, Nano Letters 18, 1104 (2018), pMID: 29369641.
  * Chen  _et al._ (2021) X. Chen, J. Zhang, C. Wen, K. Liu, Z. Zhu, S. Qin, and X. Yuan, Carbon 173, 126 (2021).
  * Cheng  _et al._ (2021) L. Cheng, R. Alaee, A. Safari, M. Karimi, L. Zhang, and R. W. Boyd, ACS Photonics 8, 585 (2021).
  * Yang  _et al._ (2017a) K.-Y. Yang, R. Verre, J. Butet, C. Yan, T. J. Antosiewicz, M. Käll, and O. J. F. Martin, Nano Letters 17, 5258 (2017a), pMID: 28829601\.
  * Xu  _et al._ (2020) L. Xu, G. Saerens, M. Timofeeva, D. A. Smirnova, I. Volkovskaya, M. Lysevych, R. Camacho-Morales, M. Cai, K. Zangeneh Kamali, L. Huang, F. Karouta, H. H. Tan, C. Jagadish, A. E. Miroshnichenko, R. Grange, D. N. Neshev, and M. Rahmani, ACS Nano 14, 1379 (2020), pMID: 31877017.
  * Nauman  _et al._ (2021) M. Nauman, J. Yan, D. de Ceglia, M. Rahmani, K. Zangeneh Kamali, C. De Angelis, A. E. Miroshnichenko, Y. Lu, and D. N. Neshev, Nature Communications 12, 1 (2021).
  * Tymchenko  _et al._ (2016) M. Tymchenko, J. S. Gomez-Diaz, J. Lee, N. Nookala, M. A. Belkin, and A. Alù, Phys. Rev. B 94, 214303 (2016).
  * Bar-David and Levy (2019) J. Bar-David and U. Levy, Nano Letters 19, 1044 (2019).
  * Minovich  _et al._ (2015) A. E. Minovich, A. E. Miroshnichenko, A. Y. Bykov, T. V. Murzina, D. N. Neshev, and Y. S. Kivshar, Laser & Photonics Reviews 9, 195 (2015).
  * Li  _et al._ (2017) G. Li, S. Zhang, and T. Zentgraf, Nature Reviews Materials 2, 1 (2017).
  * Krasnok  _et al._ (2018) A. Krasnok, M. Tymchenko, and A. Alù, Materials Today 21, 8 (2018).
  * Lee  _et al._ (2014) J. Lee, M. Tymchenko, C. Argyropoulos, P.-Y. Chen, F. Lu, F. Demmerle, G. Boehm, M.-C. Amann, A. Alu, and M. A. Belkin, Nature 511, 65 (2014).
  * Kauranen and Zayats (2012) M. Kauranen and A. V. Zayats, Nature Photonics 6, 737 (2012).
  * Butet  _et al._ (2015) J. Butet, P.-F. Brevet, and O. J. F. Martin, ACS Nano 9, 10545 (2015), pMID: 26474346\.
  * Boyd (2020) R. W. Boyd,  _Nonlinear Optics_, fourth edition ed. (Academic Press, 2020).
  * Chandrasekar  _et al._ (2015) R. Chandrasekar, N. K. Emani, A. Lagutchev, V. M. Shalaev, C. Ciracì, D. R. Smith, and A. V. Kildishev, Opt. Mater. Express 5, 2682 (2015).
  * Kruk  _et al._ (2015) S. Kruk, M. Weismann, A. Y. Bykov, E. A. Mamonov, I. A. Kolmychek, T. Murzina, N. C. Panoiu, D. N. Neshev, and Y. S. Kivshar, ACS Photonics 2, 1007 (2015).
  * Smirnova and Kivshar (2016) D. Smirnova and Y. S. Kivshar, Optica 3, 1241 (2016).
  * Bernasconi  _et al._ (2016) G. D. Bernasconi, J. Butet, and O. J. F. Martin, J. Opt. Soc. Am. B 33, 768 (2016).
  * Butet  _et al._ (2017) J. Butet, G. D. Bernasconi, M. Petit, A. Bouhelier, C. Yan, O. J. F. Martin, B. Cluzel, and O. Demichel, ACS Photonics 4, 2923 (2017).
  * Yang  _et al._ (2017b) K.-Y. Yang, J. Butet, C. Yan, G. D. Bernasconi, and O. J. F. Martin, ACS Photonics 4, 1522 (2017b).
  * Kiselev  _et al._ (2019) A. Kiselev, G. D. Bernasconi, and O. J. F. Martin, Opt. Express 27, 38708 (2019).
  * Gupta  _et al._ (2021) T. D. Gupta, L. Martin-Monier, J. Butet, K.-Y. Yang, A. Leber, C. Dong, T. Nguyen-Dang, W. Yan, O. J. Martin, and F. Sorin, Nanophotonics 10, 3465 (2021).
  * Abir  _et al._ (2022) T. Abir, M. Tal, and T. Ellenbogen, Nano Letters 22, 2712 (2022), pMID: 35369689\.
  * Li  _et al._ (2021) G.-C. Li, D. Lei, M. Qiu, W. Jin, S. Lan, and A. V. Zayats, Nature Communications 12, 4326 (2021).
  * Asadchy  _et al._ (2018) V. S. Asadchy, A. Díaz-Rubio, and S. A. Tretyakov, Nanophotonics 7, 1069 (2018).
  * Achouri and Caloz (2021) K. Achouri and C. Caloz,  _Electromagnetic Metasurfaces: Theory and Applications_ (John Wiley & Sons, 2021).
  * Achouri  _et al._ (2022a) K. Achouri, V. Tiukuvaara, and O. J. Martin, arXiv preprint arXiv:2208.12504 (2022a).
  * Overvig and Alù (2022) A. Overvig and A. Alù, Laser & Photonics Reviews 16, 2100633 (2022).
  * Poutrina and Urbas (2016) E. Poutrina and A. Urbas, Scientific reports 6, 1 (2016).
  * Mobini  _et al._ (2021) E. Mobini, R. Alaee, R. W. Boyd, and K. Dolgaleva, ACS Photonics 8, 3234 (2021).
  * Jin and Argyropoulos (2020) B. Jin and C. Argyropoulos, Phys. Rev. Applied 13, 054056 (2020).
  * Kim (2021) K.-H. Kim, Plasmonics 16, 77 (2021).
  * Liu  _et al._ (2021) W. Liu, L. Huang, J. Ding, C. Xie, Y. Luo, and W. Hong, Nanomaterials 11 (2021), 10.3390/nano11092410.
  * Abasahl  _et al._ (2021) B. Abasahl, C. Santschi, T. V. Raziman, and O. J. F. Martin, Nanotechnology 32, 475202 (2021).
  * Wang  _et al._ (2019) J. Wang, J. Butet, G. D. Bernasconi, A.-L. Baudrion, G. Lévêque, A. Horrer, A. Horneber, O. J. F. Martin, A. J. Meixner, M. Fleischer, P.-M. Adam, and D. Zhang, Nanoscale 11, 23475 (2019).
  * Wang  _et al._ (2021) J. Wang, A.-L. Baudrion, J. Béal, A. Horneber, F. Tang, J. Butet, O. J. F. Martin, A. J. Meixner, P.-M. Adam, and D. Zhang, The Journal of Chemical Physics 154, 074701 (2021).
  * Kruk  _et al._ (2022) S. S. Kruk, L. Wang, B. Sain, Z. Dong, J. Yang, T. Zentgraf, and Y. Kivshar, Nature Photonics 16, 561 (2022).
  * Czaplicki  _et al._ (2015) R. Czaplicki, J. Mäkitalo, R. Siikanen, H. Husu, J. Lehtolahti, M. Kuittinen, and M. Kauranen, Nano Letters 15, 530 (2015), pMID: 25521745.
  * Kottmann and Martin (2001) J. P. Kottmann and O. J. F. Martin, Opt. Lett. 26, 1096 (2001).
  * Gallinet  _et al._ (2010) B. Gallinet, A. M. Kern, and O. J. F. Martin, J. Opt. Soc. Am. A 27, 2261 (2010).
  * Butet  _et al._ (2013) J. Butet, B. Gallinet, K. Thyagarajan, and O. J. F. Martin, J. Opt. Soc. Am. B 30, 2970 (2013).
  * McPeak  _et al._ (2015) K. M. McPeak, S. V. Jayanti, S. J. P. Kress, S. Meyer, S. Iotti, A. Rossinelli, and D. J. Norris, ACS Photonics 2, 326 (2015).
  * Castro-Lopez  _et al._ (2011) M. Castro-Lopez, D. Brinks, R. Sapienza, and N. F. van Hulst, Nano Letters 11, 4674 (2011), pMID: 21970569\.
  * Knight  _et al._ (2014) M. W. Knight, N. S. King, L. Liu, H. O. Everitt, P. Nordlander, and N. J. Halas, ACS Nano 8, 834 (2014), pMID: 24274662.
  * Gérard and Gray (2014) D. Gérard and S. K. Gray, Journal of Physics D: Applied Physics 48, 184001 (2014).
  * Thyagarajan  _et al._ (2016) K. Thyagarajan, C. Santschi, P. Langlet, and O. J. F. Martin, Advanced Optical Materials 4, 871 (2016).
  * Krause  _et al._ (2004) D. Krause, C. W. Teplin, and C. T. Rogers, Journal of Applied Physics 96, 3626 (2004).
  * Sounas and Alù (2014) D. L. Sounas and A. Alù, Opt. Lett. 39, 4053 (2014).
  * Achouri  _et al._ (2017) K. Achouri, Y. Vahabzadeh, and C. Caloz, Opt. Express 25, 19013 (2017).
  * Achouri  _et al._ (2018) K. Achouri, G. D. Bernasconi, J. Butet, and O. J. F. Martin, IEEE Transactions on Antennas and Propagation 66, 6061 (2018).
  * Achouri  _et al._ (2022b) K. Achouri, A. Kiselev, and O. J. F. Martin, New Journal of Physics 24, 025006 (2022b).
  * Ray  _et al._ (2020) D. Ray, T. V. Raziman, C. Santschi, D. Etezadi, H. Altug, and O. J. F. Martin, Nano Letters 20, 8752 (2020), pMID: 33206533.
  * Beversluis  _et al._ (2003) M. R. Beversluis, A. Bouhelier, and L. Novotny, Phys. Rev. B 68, 115433 (2003).
  * Mühlschlegel  _et al._ (2005) P. Mühlschlegel, H.-J. Eisler, O. J. F. Martin, B. Hecht, and D. W. Pohl, Science 308, 1607 (2005).
  * Biagioni  _et al._ (2009) P. Biagioni, M. Celebrano, M. Savoini, G. Grancini, D. Brida, S. Mátéfi-Tempfli, M. Mátéfi-Tempfli, L. Duò, B. Hecht, G. Cerullo, and M. Finazzi, Phys. Rev. B 80, 045411 (2009).
  * Toyoda  _et al._ (2021) S. Toyoda, M. Fiebig, T. hisa Arima, Y. Tokura, and N. Ogawa, Science Advances 7, eabe2793 (2021).
  * Mund  _et al._ (2021) J. Mund, D. R. Yakovlev, A. N. Poddubny, R. M. Dubrovin, M. Bayer, and R. V. Pisarev, Phys. Rev. B 103, L180410 (2021).
  * Trzeciecki and Hübner (2000) M. Trzeciecki and W. Hübner, Phys. Rev. B 62, 13888 (2000).
  * Sounas and Alù (2017) D. L. Sounas and A. Alù, Phys. Rev. Lett. 118, 154302 (2017).
  * Willner  _et al._ (2014) A. E. Willner, S. Khaleghi, M. R. Chitgarha, and O. F. Yilmaz, J. Lightwave Technol. 32, 660 (2014).

## Supporting Information

Figure S1: Control over linear and SH transmission via $L_{x}$ and $L_{y}$
geometrical parameters at excitation wavelength
$\lambda_{0}=$800\text{\,}\mathrm{nm}$$. Other parameters are fixed:
$D=t_{\ce{Ag}}=t_{\ce{Al}}=$50\text{\,}\mathrm{nm}$$ and
$L_{\textrm{s}}=$25\text{\,}\mathrm{n}\mathrm{m}$$. (a) Linear and (b) SH
transmission upon forward excitation; (c) linear and (d) SH transmission upon
backward excitation. (e) Forward/backward-excitation SH transmission
extinction ratio.  Figure S2: Dependence of the linear and SH transmission and
reflection upon forward excitaiton (FT and FR) on the geometrical parameters
(a) Linear and (b) SH. Other geometrical parameters are fixed:
$L_{x}=$$135\text{\,}\mathrm{nm}$, $L_{y}=$$195\text{\,}\mathrm{nm}$ and
$D=t_{\ce{Ag}}=t_{\ce{Al}}=$50\text{\,}\mathrm{nm}$$.  Figure S3: Electric
field components calculated for the reflected and transmitted SH waves. In
both cases, (a) forward excitation (FE) and (b) backward excitation (BE) the
$E_{y}^{2\omega}$ component (that is orthogonal to the excitation field
$E_{x}^{\omega}$) is dominant.  Figure S4: Dielectric permittivity of Al (blue
lines) and Ag (gray lines) used in the simulations: real (bottom panel) and
imaginary (top panel) parts of the interpolated experimental data from ref.
McPeak  _et al._ (2015).  Figure S5: Isolated meta-atom linear scattering and
absorption analysis. (a) Scattering (red solid lines), absorption (red dotted
lines) and extinction (black-solid lines) cross-sections upon forward
excitation (upward triangles) and backward excitation (downward triangles);
(b) absoption in Ag (gray solid lines) and Al (blue dotted lines) domains upon
forward excitation (upward triangles) and backward excitation (downward
triangles). Pseud-color images of electric field distribution upon forward (c)
and backward (d) excitations. Normalized magnitude of the electric field on
the logarithmic scale is plotted.  Figure S6: Isolated meta-atom multipole
analysis. Vector spherical harmonic decomposition of (a) linear scattering
upon forward excitation; (b) linear scattering upon backward excitation; (c)
SHG scattering upon forward excitation; (d) SHG scattering upon backward
excitation.


# Betatron frequency and the Poincaré rotation number

Sergei Nagaitsev  Fermilab, Batavia, IL 60510, USA  The Enrico Fermi
Institute, The University of Chicago, Chicago, IL 60637, USA  Timofey Zolkin
Fermilab, Batavia, IL 60510, USA

###### Abstract

Symplectic maps are routinely used to describe single-particle dynamics in
circular accelerators. In the case of a linear accelerator map, the rotation
number (the betatron frequency) can be easily calculated from the map itself.
In the case of a nonlinear map, the rotation number is normally obtained
numerically, by iterating the map for given initial conditions, or through a
normal form analysis, a type of a perturbation theory for maps. Integrable
maps, a subclass of symplectic maps, allow for an analytic evaluation of their
rotation numbers. In this paper we propose an analytic expression to determine
the rotation number for integrable symplectic maps of the plane and present
several examples, relevant to accelerators. These new results can be used to
analyze the topology of the accelerator Hamiltonians as well as to serve as
the starting point for a perturbation theory for maps.

##  I Introduction

The first mention of the betatron frequency was in the 1941 pioneering work by
Kerst and Serber Kerst and Serber (1941), where they defined it as the
fractional number of particle oscillations around the orbit per one revolution
period in a betatron (a type of an induction accelerator). Later, the theory
of the alternating-gradient (AG) synchrotron Courant and Snyder (1958)
demonstrated the existence of an integral of motion (the so-called Courant-
Snyder invariant) for particles in an AG synchrotron and established a
powerful connection between the modern AG focusing systems and linear
symplectic maps, thus connecting the betatron frequency and the Poincaré
rotation number Poincaré (1885).

In modern accelerators (for example, in the LHC) particles are stored for up
to $10^{9}$ revolutions and understanding their dynamics is crucially
important for maintaining long-term particle stability Todesco (1998);
Papaphilippou (2014). One important parameter of particle dynamics in an
accelerator is the betatron frequency and its dependence on a particle’s
amplitude. It turns out that the accelerator focusing systems conserve the
Courant-Snyder invariant only approximately and there is a need to analyze the
conditions for stable particle dynamics. Over the recent years, several
methods were developed to analyze the particle motion in accelerator systems,
using either numeric tools, like the Frequency Map Analysis Dumas and Laskar
(1993), or the Normal Form Analysis Bazzani  _et al._ (1994); Turchetti and
Panichi (2019), a type of a perturbation theory, which uses a linear map and a
Courant-Snyder invariant as a starting point.

At the same time, there has been continuous interest, starting with E.
McMillan McMillan (1971), in making the accelerator maps nonlinear, yet
integrable Chow and Cary (1994); Danilov (2008); Danilov and Nagaitsev (2010,
2014). However, there does not exist an analytic method to calculate the
betatron frequency (the Poincaré rotation number) for nonlinear symplectic
integrable maps. This present paper is set to remedy this deficiency.

##  II Betatron frequency

For a one degree-of-freedom time-independent system, the Hamiltonian function,
$\mathrm{H}[p,q;t]=\mathrm{E}$, is the integral of the motion. If the motion
is bounded, it is also periodic, and the period of oscillations can be
determined by integrating

$T(\mathrm{E})=\oint\left(\frac{\partial\mathrm{H}}{\partial
p}\right)^{-1}\,\mathrm{d}q,$ (1)

where $p=p(\mathrm{E},q)$ Lichtenberg and Lieberman (1992). The oscillation
period and its dependence on initial conditions is one of the key properties
of the periodic motion.

Let us now consider a symplectic map of the plane (corresponding to a one-turn
map of an accelerator), $\mathrm{M}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$,

$(q^{\prime},p^{\prime})=\mathrm{M}\,(q,p),$

where the prime symbols (′) indicate the transformed phase space coordinates.
Suppose that the sequence, generated by a repeated application of the map,

$(q_{0},p_{0})\rightarrow(q_{1},p_{1})\rightarrow(q_{2},p_{2})\rightarrow(q_{3},p_{3})\rightarrow\ldots$

belongs to a closed invariant curve. We do not describe how this map is
obtained (see, for example, Dragt (2013)) but let us suppose that we know the
mapping equations. Let $R_{n}$ be the rotation angle in the phase space
$(q,p)$ around a stable fixed point between two consecutive iterations
$(q_{n},p_{n})$ and $(q_{n+1},p_{n+1})$. Then, the limit, when it exists,

$\nu=\lim_{N\to\infty}\frac{1}{2\,\pi\,N}\,\sum_{n=0}^{N}R_{n}$ (2)

is called the rotation number (the betatron frequency of the one-turn map) for
that particular orbit of the map $\mathrm{M}$ Dilão and Alves-Pires (1996).
Unlike Eq. (1), which allows to express the oscillation period analytically,
Eq. (2) can be only evaluated numerically for each orbit. Let us now suppose
that there exists a non-constant real-valued continuous function
$\mathcal{K}(q,p)$, which is invariant under $\mathrm{M}$. The function
$\mathcal{K}(q,p)$ is called integral and the map is called integrable. In
this paper, we are describing the case, for which the level sets
$\mathcal{K}=\mathrm{const}$ are compact closed curves (or sets of points) and
for which the identity

$\mathcal{K}(q^{\prime},p^{\prime})=\mathcal{K}(q,p)$ (3)

holds for all $(q,p)$. There are many examples of integrable maps, including
the famous McMillan map McMillan (1971), described below. The dynamics is in
many ways similar to that of a continuous system, however, Eq. (1) is not
directly applicable since the integral $\mathcal{K}(q,p)$ is not the
Hamiltonian function. Below, we will present an expression (the Danilov
theorem) to obtain the rotation number from $\mathcal{K}(q,p)$ for an
integrable map, $\mathrm{M}$.

The Arnold-Liouville theorem for integrable maps Arnold and Avez (1968);
Veselov (1991); Meiss (1992) states that (1) the action-angle variables exist
and (2) in these variables, consecutive iterations of integrable map
$\mathrm{M}$ lie on nested circles of radius $J$ and that the map can be
written in the form of a twist map,

$\begin{bmatrix}J_{n+1}\\\\[5.69046pt]
\theta_{n+1}\end{bmatrix}=\begin{bmatrix}J_{n}\\\\[5.69046pt]
\theta_{n}+2\,\pi\,\nu(J)\mod 2\,\pi\end{bmatrix},$ (4)

where $|\nu(J)|\leq 0.5$ is the rotation number, $\theta$ is the angle
variable and $J$ is the action variable, defined by the map $\mathrm{M}$ as

$J=\frac{1}{2\,\pi}\,\oint p\,\mathrm{d}q.$ (5)

Thus, in this paper, we would like to consider the following question: how
does one determine the rotation number, $\nu(\mathcal{K})$, from the known
integral, $\mathcal{K}(q,p)$, and the known integrable map, $\mathrm{M}$? In
addition, in the ”Examples” section we propose how to use this theorem when
only an approximate invariant is known.

##  III Danilov theorem

######  Theorem 1 (Danilov theorem Danilov (deceased)).

Suppose a symplectic map of the plane,

$(q^{\prime},p^{\prime})=\mathrm{M}\,(q,p),$

is integrable with the invariant (integral) $\mathcal{K}(q,p)$, then its
Poincaré rotation number is

$\nu(\mathcal{K})=\int_{q}^{q^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q\Bigg{/}\oint\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q,$ (6)

where the integrals are evaluated along the invariant curve,
$\mathcal{K}(q,p)$.

Figure 1:  Constant level sets of the integral
$\mathcal{K}(q,p)=\mathrm{const}$ (left). A particular curve representing a
level set of $\mathcal{K}$ and several iterates of the map $\mathrm{M}$
(center). A three-dimensional phase space, $(q,p)$ \+ time, of the system (7)
(right). Dark gray planes $t=0,\tau,2\tau,\ldots$ represent stroboscopic
Poincaré section of the continuous flow of the system (red curve) which is
identical to map $\mathrm{M}$.

###### Proof.

Consider the following system of differential equations:

$\frac{\mathrm{d}Q}{\mathrm{d}t}=\frac{\partial\mathcal{K}(Q,P)}{\partial
P},\qquad\qquad\frac{\mathrm{d}P}{\mathrm{d}t}=-\frac{\partial\mathcal{K}(Q,P)}{\partial
Q}.$ (7)

We notice that $\mathcal{K}(Q,P)$ does not change along a solution of the
system, because it is an integral of the motion, meaning

$\frac{\mathrm{d}\mathcal{K}}{\mathrm{d}t}=\frac{\partial\mathcal{K}}{\partial
Q}\frac{\mathrm{d}Q}{\mathrm{d}t}+\frac{\partial\mathcal{K}}{\partial
P}\frac{\mathrm{d}P}{\mathrm{d}t}=0$ (8)

for any solution $Q(t)$ and $P(t)$. Let $q(t)$ and $p(t)$ be the solutions of
the system (7) with the following initial conditions $q(0)=q_{0}$ and
$p(0)=p_{0}$. Define a new map, $\widetilde{\mathrm{M}}(q,p)$ (see Fig. 1)

$(q^{\prime},p^{\prime})=\widetilde{\mathrm{M}}(q,p)=\left(q(\tau),p(\tau)\right)$
(9)

where $\tau$ is a discrete time step. For a given $\mathcal{K}$, which is an
integral of both $\mathrm{M}$ and $\widetilde{\mathrm{M}}$, one can always
select $\tau(\mathcal{K})$ such that the maps $\mathrm{M}(q,p)$ and
$\widetilde{\mathrm{M}}(q,p)$ are identical. This follows from the Arnold-
Liouville theorem. Since $\mathcal{K}(q,p)$ is compact and closed, the
functions $q(t)$ and $p(t)$ are periodic with a period $T(\mathcal{K})$. By
its definition,

$\tau=\nu(\mathcal{K})\,T(\mathcal{K}).$ (10)

Let us now calculate $\nu(\mathcal{K})$:

$\begin{array}[]{l}\displaystyle\nu(\mathcal{K})\equiv\frac{\tau}{T}=\frac{\int_{q}^{q^{\prime}}\,\mathrm{d}t}{\oint\,\mathrm{d}t}=\frac{\int_{q}^{q^{\prime}}\left(\frac{\mathrm{d}q}{\mathrm{d}t}\right)^{-1}\,\mathrm{d}q}{\oint\left(\frac{\mathrm{d}q}{\mathrm{d}t}\right)^{-1}\,\mathrm{d}q}\\\\[18.49411pt]
\displaystyle\qquad\,\,=\frac{\int_{q}^{q^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q}{\oint\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q}.\end{array}$ (11)

Q.E.D. ∎

Figure 2:  The partial action is defined as a sector area (blue) for one map
iteration, divided by ${2\,\pi}$ (a.). Convenient choices of the partial
action for mappings in McMillan form: an area under the curve in II (blue) or
IV (green) quadrants (b.), and, areas for initial conditions in a form of
$(q_{0},q_{0})$ (c.).

######  Corollary 1.1.

$\nu(\mathcal{K})=\frac{\mathrm{d}J^{\prime}}{\mathrm{d}J},$ (12)

where

$J^{\prime}(\mathcal{K})=\frac{1}{2\,\pi}\,\int_{q}^{q^{\prime}}p(\mathcal{K},q)\,\mathrm{d}q.$
(13)

is the partial action calculated as a sector integral (see Fig. 2) around the
stable fixed point.

###### Proof.

First, we will consider the denominator in Eq. (11):

$\frac{1}{2\,\pi}\,\oint\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q=\frac{1}{2\,\pi}\,\frac{\mathrm{d}}{\mathrm{d}\mathcal{K}}\,\oint
p\,\mathrm{d}q=\frac{\mathrm{d}J}{\mathrm{d}\mathcal{K}}.$ (14)

Second, we will evaluate the numerator. Using the equations of motion in Eq.
(7), we notice that

$\int_{q}^{q^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q=-\int_{p}^{p^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
q}\right)^{-1}\,\mathrm{d}p.$ (15)

Now, we will utilize the Leibniz integral rule together with Eq. (15) to
obtain

$\begin{array}[]{l}\displaystyle\frac{1}{2\,\pi}\,\int_{q}^{q^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p}\right)^{-1}\,\mathrm{d}q=\frac{1}{2\,\pi}\times\\\\[18.49411pt]
\displaystyle\,\,\,\,\,\,\frac{\mathrm{d}}{\mathrm{d}\mathcal{K}}\left(\frac{q\,p-q^{\prime}\,p^{\prime}}{2}+\int_{q}^{q^{\prime}}p\,\mathrm{d}q\right)=\frac{\mathrm{d}J^{\prime}}{\mathrm{d}\mathcal{K}}.\end{array}$
(16)

Finally, by combining Eqs. (14) and (16) we obtain the Eq. (12). ∎

######  Corollary 1.2.

For a linear map ($\nu=\mathrm{const}$),

$\nu=J^{\prime}/J.$ (17)

###### Proof.

Since $\nu=\mathrm{const}$, the Hamiltonian function is
$\mathrm{H}(J)=\nu\,J$. Using Eq. (12), we obtain Eq. (17). ∎

######  Corollary 1.3.

The Hamiltonian function corresponding to the map $M$ is

$\mathrm{H}(\mathcal{K})=J^{\prime}(\mathcal{K}).$ (18)

###### Proof.

Since $\nu=\mathrm{d}\mathrm{H}/\mathrm{d}J$, one can use Eq. (12) to obtain
$\mathrm{H}=J^{\prime}+\mathrm{const}$. ∎

######  Corollary 1.4.

$\nu(\mathcal{K})=\int_{p}^{p^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
q}\right)^{-1}\,\mathrm{d}p\Bigg{/}\oint\left(\frac{\partial\mathcal{K}}{\partial
q}\right)^{-1}\,\mathrm{d}p,$ (19)

where the integrals are evaluated along the invariant curve,
$\mathcal{K}(q,p)$.

###### Proof.

Because of the $p\leftrightarrow-q$ symmetry in Eqs. (7), the proof is similar
to Eq. (11). ∎

In order to generalize the Danilov theorem to higher-dimensional integrable
maps, one has to know the variables, where such a map is separated into maps
for each degree of freedom. Below we will consider an example of a 4D map,
which is separable in polar coordinates with two integrals of motion.

##  IV Examples

In order to employ this theorem in practice, one would need to recall that
with $p=p(\mathcal{K},q)$, the integrand in Eq. (6) is

$\left(\frac{\partial\mathcal{K}}{\partial p}\right)^{-1}=\frac{\partial
p(\mathcal{K},q)}{\partial\mathcal{K}}.$ (20)

Also, the lower limit of the integral can be chosen to be any convenient value
of $q$, for example 0, as long it belongs to a given level set,
$\mathcal{K}(q,p)$. Finally, the upper limit of the integral, $q^{\prime}$, is
obtained from the selected $q$ and $p=p(\mathcal{K},q)$ by iterating the map,
$\mathrm{M}(q,p)$. It is clear that not all functions $\mathcal{K}(q,p)$ can
be inverted analytically to obtain $p=p(\mathcal{K},q)$. This drawback of this
method can be overcome by numeric evaluations (see Appendix B).

For maps in a special (McMillan) form McMillan (1971),

$\begin{bmatrix}q^{\prime}\\\ p^{\prime}\end{bmatrix}=\begin{bmatrix}p\\\
-q+f(p)\end{bmatrix},$ (21)

the convenient choices for integration limits in Eq. 6 are $(q,p)=(q_{0},0)$
and $(q^{\prime},p^{\prime})=(0,-q_{0}+f(0))$, Fig. 2.b, and $(q,p)=(a,a)$ and
$(q^{\prime},p^{\prime})=(a,-a+f(a))$, Fig. 2.c. Finally, for twist maps, Eq.
(4), the Danilov theorem Eq. (6) gives $\nu$, as expected.

Let us now consider several non-trivial examples. Linear maps are presented in
Appendix A.

###  IV.1 McMillan map

As our first example, we will consider the so-called McMillan map McMillan
(1971),

$\begin{bmatrix}q^{\prime}\\\ p^{\prime}\end{bmatrix}=\begin{bmatrix}p\\\
-q+a\,p/(b\,p^{2}+1)\end{bmatrix}.$ (22)

This map has been considered in detailed in Ref. Iatrou and Roberts (2002). To
illustrate the Danilov theorem, we will limit ourselves to a case with $b>0$
and $|a|<2$, which corresponds to stable motion at small amplitudes. Mapping
(22) has the following integral:

$\mathcal{K}(q,p)=b\,q^{2}p^{2}+q^{2}+p^{2}-a\,q\,p,$ (23)

which is non-negative for the chosen parameters.

We first notice that for small amplitudes, $b\,p^{2}\ll 1$, this map can be
approximated as

$\begin{bmatrix}q^{\prime}\\\
p^{\prime}\end{bmatrix}\approx\begin{bmatrix}p\\\
-q+a\,p-a\,b\,p^{3}+a\,b^{2}\,p^{5}-...\end{bmatrix},$ (24)

and its zero-amplitude rotation number is Courant and Snyder (1958)

$\nu(0)=\frac{1}{2\,\pi}\,\arccos\frac{a}{2}.$ (25)

At large amplitudes ($b\,p^{2}\gg 1$), the rotation number becomes $0.25$. We
will now evaluate the rotation number analytically, using Eq. (6): Let us
define a parameter,

$w(\mathcal{K})=\frac{1}{\sqrt{2}}\sqrt{1+\frac{d(\mathcal{K})}{\sqrt{d(\mathcal{K})^{2}+4\,\mathcal{K}\,b}}},$
(26)

which spans from 0 to 1 and where $d(\mathcal{K})=a^{2}/4+\mathcal{K}\,b-1$.
Then, the rotation number can be expressed through Jacobi elliptic functions
as follows:

$\nu(\mathcal{K})=\frac{1}{4\,\mathrm{K}(w)}\mathrm{arcds}\left({\left(d(\mathcal{K})^{2}+4\,\mathcal{K}\,b\right)^{-1/4}},w\right),$
(27)

where $\mathrm{K}(w)$ is the complete elliptic integral of the first kind and
the inverse Jacobi function, $\mathrm{arcds}(x,w)$, is defined as follows

$\mathrm{arcds}(x,w)=\int_{x}^{\infty}\frac{\mathrm{d}t}{\sqrt{(t^{2}+w^{2})(t^{2}+w^{2}-1)}}.$
(28)

The rotation number, Eq. (27), has the following series expansion:

$\nu(\mathcal{K})\approx\nu(0)+\frac{3}{2\,\pi}\,\frac{b\,a}{\sqrt{(4-a^{2})^{3}}}\mathcal{K}.$
(29)

Figure 3:  The left plot contains iterations (green dots) of the McMillan map
($a=1.6$, $b=1$). Constant level sets of the invariant are shown with blue
lines. The right plot is the rotation number, Eq. (27), as a function of its
integral, $\mathcal{K}$. The inset shows the linear approximation, Eq. (29).

Figure 3 shows an example of the rotation number, for the case of $a=1.6$ and
$b=1$ ($\nu(0)\approx 0.102$), as a function of integral, $\mathcal{K}$.

The McMillan invariant (23) also allows for an analytic evaluation of the
action integral (5). We will omit the lengthy expressions, but will only
present a small-amplitude series expansion:

$J(\mathcal{K})\approx\frac{\mathcal{K}}{\sqrt{4-a^{2}}}-\frac{b\,(2+a^{2})\,\mathcal{K}^{2}}{\sqrt{(4-a^{2})^{5}}}.$
(30)

Finally, we can also present a small-amplitude series expansion of the
rotation number (27):

$\nu(J)\approx\nu(0)+\frac{3}{2\,\pi}\,\frac{b\,a}{4-a^{2}}J.$ (31)

###  IV.2 Cubic map

Figure 4:  The top row: phase-space trajectories of a cubic map, obtained by
tracking with $a=-0.85$ (left plot) and level sets of the approximate
invariant (34) (right plot), on the same scale. The red and blue lines in the
top left plot corresponds to symmetry lines $p=q$ and
$p=(a\,q+\epsilon\,q^{3})/2$ respectively. The bottom row: the left plot shows
the rotation number as a function of initial conditions in the form
$q_{0}=p_{0}$, by using the Eq. (2) (black solid line), and by using the
Danilov theorem, Eq. (12) numerically (orange dashed). The red solid line
corresponds to the rotation number obtained from the approximate invariant
(34) using the Danilov theorem as well. The right bottom plot shows the
dependence of $\nu$ as a function of action $J$, from tracking (orange dashed)
and from the approximate invariant (34) (red solid).

As our second example, we will consider a non-integrable Hénon cubic map Hénon
(1969); Dullin and Meiss (2000):

$\begin{bmatrix}q^{\prime}\\\ p^{\prime}\end{bmatrix}=\begin{bmatrix}p\\\
-q+a\,p+\epsilon\,p^{3}\end{bmatrix}.$ (32)

This map is well-known in accelerator physics as a symplectic octupole map. At
small amplitudes this map is linear and the rotation number is

$\nu\approx\frac{1}{2\,\pi}\,\arccos\left(\frac{a}{2}\right).$ (33)

At large amplitudes this map becomes chaotic and unstable. Let us propose an
approximate integral (the exact integral does not exist since it is a non-
integrable map).

$\begin{array}[]{l}\displaystyle\mathcal{K}_{\text{c}}(q,p)=p^{2}+q^{2}-a\,p\,q-\frac{\epsilon}{a}\,p^{2}q^{2}\\\\[8.5359pt]
\displaystyle\qquad+\frac{7\,\epsilon}{5\,a\,(4-a^{2})}\,\left(p^{2}+q^{2}-a\,p\,q\right)^{2}+O\left(\epsilon^{2}\right).\end{array}$
(34)

The derivation of this approximate integral goes beyond the scope of this
article and will be described in subsequent publications. For this
illustration, the reader can verify by inspection that this integral is
approximately conserved, near the origin. We will now use the Danilov theorem
to evaluate the rotation number of this map for various initial conditions
with $q_{0}=p_{0}$. Figure 4 shows the exact (numeric), Eq. (2), and the
approximate rotation number, calculated from (34) and (32) using the Danilov
theorem, Eq. (6).

A small-amplitude series expansion of the rotation number is:

$\nu(J)\approx\nu(0)-\frac{3}{2\pi}\frac{\epsilon}{4-a^{2}}J,$ (35)

which is the same as in Dullin and Meiss (2000) and similar to Eq. (31).

###  IV.3 4-D integrable map

In this section we will sketch out an example of how to use the Danilov
theorem to analyze an integrable multi-dimensional map. Consider the following
map, which can be realized in accelerators by employing the so-called electron
lens Shiltsev  _et al._ (1999, 2008); Lobach  _et al._ (2018),

$\begin{bmatrix}x^{\prime}\\\\[1.42271pt] p_{x}^{\prime}\\\\[1.42271pt]
y^{\prime}\\\\[1.42271pt]
p_{y}^{\prime}\end{bmatrix}=\begin{bmatrix}\alpha_{x}x+\beta\,p_{x}\\\
-\gamma_{x}x-\alpha_{x}\,p_{x}+\frac{a\,x^{\prime}}{b\,r^{\prime 2}+1}\\\
\alpha_{y}y+\beta\,p_{y}\\\
-\gamma_{y}y-\alpha_{y}\,p_{y}+\frac{a\,y^{\prime}}{b\,r^{\prime
2}+1}\end{bmatrix},$ (36)

where $r^{2}=x^{2}+y^{2}$, $\beta\,\gamma_{x}=1+\alpha_{x}^{2}$,
$\beta\,\gamma_{y}=1+\alpha_{y}^{2}$, with $\alpha_{x}$, $\alpha_{y}$, $a$,
$b$ and $\beta$ being some arbitrary parameters. This map has two integrals of
motion in involution (having a vanishing Poisson bracket):

$L=(\alpha_{y}-\alpha_{x})\,x\,y+\beta(x\,p_{y}-y\,p_{x})$ (37)

and

$\mathcal{K}=\Big{(}b+\frac{1}{r^{2}}\Big{)}\,T^{2}+\beta\,a\,T+r^{2}+\frac{L^{2}}{r^{2}},$
(38)

where $T=\alpha_{x}x^{2}+\alpha_{y}y^{2}+\beta\,r\,p_{r}$ and
$p_{r}=(x\,p_{x}+y\,p_{y})/r$. In order to employ the Danilov theorem, we must
rewrite the map (36) in new variables, where this map is separated into two
maps. Such variables exist by virtue of this map being integrable. We first
notice that by introducing new variables,

$\begin{matrix}\tilde{x}=x/\sqrt{\beta}\\\
\tilde{p}_{x}=x\,\alpha_{x}/\sqrt{\beta}+p_{x}\sqrt{\beta}\\\
\tilde{y}=y/\sqrt{\beta}\\\
\tilde{p}_{y}=y\,\alpha_{y}/\sqrt{\beta}+p_{y}\sqrt{\beta},\end{matrix}$ (39)

the map (36) becomes symmetric in $\tilde{x}$ and $\tilde{y}$ with
$\tilde{a}=a\sqrt{\beta}$ and $\tilde{b}=b\beta$. The resulting map is
separable in polar coordinates, $r$ and $\theta$, such that
$x=r\,\cos(\theta)$ and $y=r\,\sin(\theta)$, where we omitted the tilde
($\,\tilde{}\,$) sign for clarity. The resulting map is

$\begin{bmatrix}r^{\prime}\\\\[8.5359pt] p_{r}^{\prime}\\\\[8.5359pt]
\theta^{\prime}\\\\[8.5359pt]
p_{\theta}^{\prime}\end{bmatrix}=\begin{bmatrix}\sqrt{p_{r}^{2}+\frac{p_{\theta}^{2}}{r^{2}}}\\\\[7.11317pt]
-p_{r}\frac{r}{r^{\prime}}+\frac{a\,r^{\prime}}{b\,r^{\prime
2}+1}\\\\[7.11317pt] \theta+\arctan\frac{p_{\theta}}{r\,p_{r}}\\\\[7.11317pt]
p_{\theta}\end{bmatrix},$ (40)

where the angular momentum $p_{\theta}=x\,p_{y}-y\,p_{x}=\mathrm{const}$ is
the integral of the motion. An additional integral is

$\mathcal{K}(r,p_{r},p_{\theta})=b\,r^{2}p_{r}^{2}+r^{2}+p_{r}^{2}-a\,r\,p_{r}+\frac{p_{\theta}^{2}}{r^{2}}.$
(41)

Now we will use the Danilov theorem to obtain two unknown rotation numbers,
$\nu_{\theta}$ and $\nu_{r}$. We first notice that $\mathcal{K}$ does not
depend on $\theta$ and thus can be used to evaluate $\nu_{r}$ in Eq. (11)
directly, by treating $p_{\theta}$ as a parameter.

$\begin{array}[]{l}\displaystyle\nu_{r}(\mathcal{K},p_{\theta})=\frac{\tau}{T_{r}}=\frac{\int_{r}^{r^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p_{r}}\right)^{-1}\,\mathrm{d}r}{\oint\left(\frac{\partial\mathcal{K}}{\partial
p_{r}}\right)^{-1}\,\mathrm{d}r}\\\\[22.76228pt]
\displaystyle\qquad=\mathrm{F}\left[\arcsin\sqrt{\frac{\zeta_{3}-\zeta_{1}}{\zeta_{3}+1}},\kappa\right]\Big{/}\left(2\,\mathrm{K}\left(\kappa\right)\right),\end{array}$
(42)

where $\mathrm{K}(\kappa)$ is the complete elliptic integral of the first
kind, $\mathrm{F}(\phi,\kappa)$ is the incomplete elliptic integral of the
first kind, elliptic modulus $\kappa$ is given by

$\kappa=\sqrt{\frac{\zeta_{3}-\zeta_{2}}{\zeta_{3}-\zeta_{1}}},$

and $\zeta_{1}<0<\zeta_{2}<\zeta_{3}$ are the roots of the polynomial

$\mathcal{P}_{3}(\zeta)=-\zeta^{3}+\left[\mathcal{K}+\left(\frac{a}{2}\right)^{2}-1\right]\,\zeta^{2}+(\mathcal{K}-p_{\theta}^{2})\,\zeta-
p_{\theta}^{2}.$

In order to evaluate the angular rotation number, $\nu_{\theta}$, we first
notice that there is some uncertainty as to which integral of the motion to
employ: one can add an arbitrary function of $p_{\theta}$ to $\mathcal{K}$,
$\mathcal{K}^{\prime}=\mathcal{K}+f\left(p_{\theta}\right)$, to obtain another
integral. This new integral of motion, $\mathcal{K}^{\prime}$, gives the same
$\nu_{r}$, but modifies the angular motion by some unknown linear function of
time:

$\frac{d\theta}{dt}=\frac{d\mathcal{K}^{\prime}}{dp_{\theta}}=\frac{d\mathcal{K}}{dp_{\theta}}+f^{\prime}\left(p_{\theta}\right),$
(43)
$\theta(t)=\int\frac{d\mathcal{K}}{dp_{\theta}}dt+f^{\prime}\left(p_{\theta}\right)t.$
(44)

Fortunately, we can resolve this uncertainty by using the angular portion of
the map, Eq. (40). By its definition, the angular rotation number is

$\nu_{\theta}=\nu_{r}\frac{\Delta\theta\left(T_{r}\right)}{2\,\pi},$ (45)

where

$\Delta\theta\left(T_{r}\right)=\oint\frac{d\mathcal{K}}{dp_{\theta}}\left(\frac{\partial\mathcal{K}}{\partial{p_{r}}}\right)^{-1}dr+k\,T_{r},$
(46)

$k$ is an unknown coefficient and $T_{r}$ is the period of the radial motion,

$T_{r}=\oint\left(\frac{\partial\mathcal{K}}{\partial
p_{r}}\right)^{-1}\,\mathrm{d}r.$ (47)

To determine the coefficient $k$ we will notice from Eq. (40) that
$\Delta\theta(\tau)=\arctan\left(\frac{p_{\theta}}{r\,p_{r}}\right)$. Thus,

$k=\frac{1}{\tau}\left(\arctan\left(\frac{p_{\theta}}{r\,p_{r}}\right)-\int_{r}^{r^{\prime}}\frac{d\mathcal{K}}{dp_{\theta}}\left(\frac{\partial\mathcal{K}}{\partial{p_{r}}}\right)^{-1}dr\right)$
(48)

with

$\tau=\int_{r}^{r^{\prime}}\left(\frac{\partial\mathcal{K}}{\partial
p_{r}}\right)^{-1}\,\mathrm{d}r.$ (49)

Now, recalling that $\nu_{r}=\tau/{T_{r}}$, we finally obtain

$\begin{array}[]{l}\displaystyle\nu_{\theta}=\frac{\nu_{r}}{2\,\pi}\oint\frac{d\mathcal{K}}{dp_{\theta}}\left(\frac{\partial\mathcal{K}}{\partial{p_{r}}}\right)^{-1}dr\,+\\\\[9.95863pt]
\displaystyle\frac{1}{2\pi}\left(\arctan\left(\frac{p_{\theta}}{r\,p_{r}}\right)-\int_{r}^{r^{\prime}}\frac{d\mathcal{K}}{dp_{\theta}}\left(\frac{\partial\mathcal{K}}{\partial{p_{r}}}\right)^{-1}dr\right).\end{array}$
(50)

After some math, this expression can be rewritten as

$\begin{array}[]{l}\displaystyle\nu_{\theta}(\mathcal{K},p_{\theta})=\frac{\Delta}{2\,\pi}\,\left[\nu_{r}-\frac{\Delta^{\prime}}{\Delta}+\frac{\arctan\left(\frac{2\,p_{\theta}}{a}\,\frac{\zeta_{3}+1}{\zeta_{3}}\right)}{\Delta},\right]\end{array}$
(51)

where

$\begin{array}[]{l}\displaystyle\Delta=\frac{2\,p_{\theta}}{\zeta_{3}\,\sqrt{\zeta_{3}-\zeta_{1}}}\,\Pi\left[\kappa\,\bigg{|}\frac{\zeta_{3}-\zeta_{2}}{\zeta_{3}}\right],\\\\[9.95863pt]
\displaystyle\Delta^{\prime}=\frac{p_{\theta}}{\zeta_{3}\,\sqrt{\zeta_{3}-\zeta_{1}}}\,\Pi\left[\arcsin\sqrt{\frac{\zeta_{3}-\zeta_{1}}{\zeta_{3}+1}},\kappa\,\Bigg{|}\frac{\zeta_{3}-\zeta_{2}}{\zeta_{3}}\right],\end{array}$

and, $\Pi(\kappa\,|\alpha)$ and $\Pi(\phi,\kappa\,|\alpha)$ are the complete
and the incomplete elliptic integrals of the third kind, respectively. One can
note that for a linear 4D map ($b=0$), we have $\nu_{r}=2\,\nu_{\theta}$ for
any value of $p_{\theta}$. Fig. 5 shows an example of the radial and the
angular rotation numbers as a function of $\mathcal{K}$ for various values of
$p_{\theta}$.

Figure 5:  Radial (left) and angular (right) rotation numbers as a function of
the first integral of the map, $\mathcal{K}$, for different values of its
second integral, $p_{\theta}$ (shown with color labels). The map parameters
are $b=1$ and $a=1.6$. Note that for $p_{\theta}=0$,
$\nu_{r}=2\,\nu_{\theta}$, as expected, and equals to the frequency $\nu$ from
the one-dimensional example of Fig. 3

##  V Summary

In this paper we demonstrated a general and exact method of how to find a
Poincaré rotation number for integrable symplectic maps of a plane and its
connection to accelerator physics. It complements the discrete Arnold-
Liouville theorem for maps Veselov (1991); Arnold and Avez (1968) and permits
the analysis of dynamics for integrable systems. Eq. (18) also permits to
express the Hamilton function of a given integrable map explicitly. Several
examples were presented in our paper. These examples demonstrate that the
Danilov theorem is a powerful tool. The McMillan integrable map is a classic
example of a nonlinear integrable discrete-time system, which finds
applications in many areas of physics, including accelerators Antipov  _et
al._ (2017); Lobach  _et al._ (2018). It is a typical member of a wide class
of area-preserving transformations called a twist map Meiss (1992). For non-
integrable maps, which are also very common in accelerator science, this new
theorem could allow for an approximate evaluation of rotation numbers,
provided there exists an approximate integral of motion, like Eq. (34).

##  VI Acknowledgments

The authors would like to thank Jeffrey Holmes and Stanislav Baturin for
carefully reading this manuscript and for their helpful comments. This
research is supported by Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy and by the University of
Chicago.

##  Appendix A Linear maps

In this appendix we will consider two examples of linear maps and we will use
Eq. (6) for one and Eq. (17) for the second one.

###  A.1 Linear accelerator map

Consider a linear symplectic map,

$\begin{bmatrix}q^{\prime}\\\ p^{\prime}\end{bmatrix}=\begin{bmatrix}a&b\\\
c&d\end{bmatrix}\begin{bmatrix}q\\\ p\end{bmatrix},$ (52)

with $a\,d-b\,c=1$ and $|a+d|\leq 2$. This map is very common in accelerator
physics and has been described in Courant and Snyder (1958). The rotation
number (the betatron frequency) for this map is well known:

$\nu=\frac{1}{2\,\pi}\,\arccos\frac{a+d}{2}.$ (53)

To obtain this equation using the Danilov theorem, we will recall that this
map has the following Courant-Snyder integral (invariant):

$\mathcal{K}=c\,q^{2}+(d-a)\,q\,p-b\,p^{2}.$ (54)

Let us assume that $c>0$, then $b\leq 0$ and $\mathcal{K}(q,p)\geq 0$ for any
$q$ and $p$. From this, we obtain

$\left(\frac{\partial\mathcal{K}}{\partial p}\right)^{-1}=\frac{\partial
p}{\partial\mathcal{K}}=\frac{\pm
1}{\sqrt{\left[(a+d)^{2}-4\right]\,q^{2}-4\,b\,\mathcal{K}}}.$ (55)

We will use

$(q,p)=(\sqrt{\mathcal{K}/c},0)$ (56)

and

$(q^{\prime},p^{\prime})=(b\,\sqrt{\mathcal{K}/c},\sqrt{\mathcal{K}\,c})$ (57)

After a straightforward evaluation of integrals in Eq. (6), we obtain:

$\nu=\frac{1}{2\,\pi}\,\arccos\frac{a+d}{2},$ (58)

same is in Eq. (53).

###  A.2 Brown map

As a second example we will consider the Brown map Brown (1983, 1993),
$\mathrm{M}_{\mathrm{B}}$,

$\begin{bmatrix}q^{\prime}\\\ p^{\prime}\end{bmatrix}=\begin{bmatrix}p\\\
-q+|p|\end{bmatrix},$ (59)

which has the following integral,

$\displaystyle\mathcal{K}(q,p)$ $\displaystyle=$
$\displaystyle\frac{1}{8}\Bigg{(}q+p+\big{|}q-|p|\big{|}+\big{|}p-|q|\big{|}+$
(60) $\displaystyle
2\,\Big{|}q-\big{|}p-|q|\big{|}\Big{|}+2\,\Big{|}p-\big{|}q-|p|\big{|}\Big{|}+$
$\displaystyle\bigg{|}q-|p|+\Big{|}p-\big{|}q-|p|\big{|}\Big{|}\bigg{|}+$
$\displaystyle\bigg{|}p-|q|+\Big{|}q-\big{|}p-|q|\big{|}\Big{|}\bigg{|}\Bigg{)}.$

The map has only one stable fixed point, located at the origin, with
$\mathcal{K}=0$. Constant level sets of $\mathcal{K}>0$ are polygons,
geometrically similar to each other, with 9 sides, labeled by Roman numerals,
see Fig. 6.a. All orbits belonging to these levels are periodic with

$\mathrm{M}_{\mathrm{B}}^{9}(q,p)=(q,p),$ (61)

and in fact, they are permutation 9-cycles such that

$\displaystyle\ldots\rightarrow\text{I}$ $\displaystyle\rightarrow$
$\displaystyle\text{III}\rightarrow\text{V}\rightarrow\text{VII}\rightarrow\text{IX}\rightarrow$
$\displaystyle\rightarrow$
$\displaystyle\text{II}\rightarrow\text{IV}\rightarrow\text{VI}\rightarrow\text{VIII}\rightarrow\text{I}\rightarrow\ldots.$

Figure 6: Brown map. (a.) Constant level sets of the invariant,
$\mathcal{K}(q,p)=\mathrm{const}$ (black solid polygons). Dashed black line
$p=q$ and blue line $p=\frac{1}{2}\,|q|$ illustrate two reflection symmetries
of the invariant polygons. Line segments are labeled with Roman numerals.
Green points are an example of a 9-cycle orbit, where the Arabic numerals show
the iteration number. (b.) An example of possible contour of integration for
the numerator and denominator in Danilov theorem.

Since it is a linear map ($\nu=\mathrm{const}$ for all orbits), we will use
Eq. (17) to determine its rotation number. It is obvious from Fig. 6.b that
$J=4.5\,\alpha$, while $J^{\prime}=1\,\alpha$, where $\alpha$ is some
arbitrary scale parameter, resulting in $\nu=\frac{2}{9}$.

##  Appendix B Numerical procedure for Danilov theorem

In this appendix we will consider two numerical procedures, which can be
employed in order to use Eq. (12) for mappings in McMillan form when only the
mapping equation is known or when we have an approximate (or an exact)
invariant of the motion but we can not compute action integrals analytically.

We will start with the case when we have only the mapping equations. As a
first step we will rewrite the map in polar coordinates

$q=r\,\cos\phi,\qquad\qquad p=r\,\sin\phi.$

Then we will iterate for various initial conditions $q_{\text{ini}}^{(k)}$,
let say in a form of $q_{\text{ini}}^{(k)}=q_{0}^{(k)}=p_{0}^{(k)}$, so that
we have a collection of points in a form

$(r_{0}^{(k)},\phi_{0}^{(k)}),(r_{1}^{(k)},\phi_{1}^{(k)}),(r_{2}^{(k)},\phi_{2}^{(k)}),\ldots,(r_{n}^{(k)},\phi_{n}^{(k)}).$

We can then sort each orbit such that

$\widetilde{\phi}_{0}^{(k)}<\widetilde{\phi}_{1}^{(k)}<\widetilde{\phi}_{2}^{(k)}<\ldots<\widetilde{\phi}_{n}^{(k)},\vspace{0.1cm}$

where $(\widetilde{r}_{i}^{(k)},\widetilde{\phi}_{i}^{(k)})$ are the points of
a new sorted $k$-th orbit. Now, for each orbit we can compute the action and
the partial action numerically as

$J^{(k)}=\frac{1}{2\,\pi}\,\sum_{i=0}^{n}\frac{\left(\widetilde{r}_{i}^{(k)}\right)^{2}}{2}\,\left[\widetilde{\phi}_{i}^{(k)}-\widetilde{\phi}_{i-1}^{(k)}\right]$
(62)

and

$J^{\prime(k)}=\frac{1}{2\,\pi}\,\sum_{\pi/2<\widetilde{\phi}_{i}^{(k)}<\pi}\frac{\left(\widetilde{r}_{i}^{(k)}\right)^{2}}{2}\,\left[\widetilde{\phi}_{i}^{(k)}-\widetilde{\phi}_{i-1}^{(k)}\right]$
(63)

respectively. Finally, using the Danilov theorem, we can find the rotation
number as a numerical derivative

$\nu^{(k)}=\frac{J^{\prime(k+1)}-J^{\prime(k)}}{J^{(k+1)}-J^{(k)}}.$ (64)

If one would like to apply the Danilov theorem directly to an approximate or
exact invariant of motion, we can proceed in a similar manner. First, we
rewrite the invariant of motion in polar coordinates,
$\mathcal{K}_{\text{approx}}(r,\phi)$. Then, for different values
$\mathcal{K}_{\text{approx}}^{(k)}$ we will numerically solve $n$ equations

$\mathcal{K}_{\text{approx}}(r,\phi^{(k)}_{i})=\mathcal{K}_{\text{approx}}^{(k)}$
(65)

with $\phi^{(k)}_{i}=2\,\pi\,i/n$ and $i=0,1,\ldots,n-1$. Denoting the
smallest positive root of equation above as $r^{(k)}_{i}$, we can find action
and partial actions as

$J^{(k)}=\frac{1}{2\,\pi}\,\sum_{i=0}^{n-1}\frac{\left(r_{i}^{(k)}\right)^{2}}{2}\,\left[\phi_{i}^{(k)}-\phi_{i-1}^{(k)}\right]$
(66)

and

$J^{\prime(k)}=\frac{1}{2\,\pi}\,\sum_{\pi/2<\phi_{i}^{(k)}<\pi}\frac{\left(r_{i}^{(k)}\right)^{2}}{2}\,\left[\phi_{i}^{(k)}-\phi_{i-1}^{(k)}\right],$
(67)

along with the rotation number

$\nu^{(k)}=\frac{J^{\prime(k+1)}-J^{\prime(k)}}{J^{(k+1)}-J^{(k)}}.$ (68)

## References

  * Kerst and Serber (1941) D. W. Kerst and R. Serber, Phys. Rev. 60, 53 (1941).
  * Courant and Snyder (1958) E. D. Courant and H. S. Snyder, Annals of physics 3, 1 (1958).
  * Poincaré (1885) H. Poincaré, Journal de mathématiques pures et appliquées 4e série 1, 167 (1885).
  * Todesco (1998) E. Todesco, in  _Analysis and Modelling of Discrete Dynamical Systems_ , edited by D. Benest and C. Froeschle (Taylor & Francis, 1998).
  * Papaphilippou (2014) Y. Papaphilippou, Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 024412 (2014).
  * Dumas and Laskar (1993) H. S. Dumas and J. Laskar, Phys. Rev. Lett. 70, 2975 (1993).
  * Bazzani  _et al._ (1994) A. Bazzani, E. Todesco, G. Turchetti, and G. Servizi, Reports- CERN (1994).
  * Turchetti and Panichi (2019) G. Turchetti and F. Panichi, “Birkhoff normal forms and stability indicators for betatronic motion,” in  _Nonlinear Dynamics and Collective Effects in Particle Beam Physics_ (World Scientific, 2019) pp. 47–69.
  * McMillan (1971) E. M. McMillan, in  _Topics in modern physics, a tribute to E.V. Condon_ , edited by E. Brittin and H. Odabasi (Colorado Assoc. Univ. Press, Boulder, CO, 1971) pp. 219–244.
  * Chow and Cary (1994) C. C. Chow and J. R. Cary, Phys. Rev. Lett. 72, 1196 (1994).
  * Danilov (2008) V. Danilov, Phys. Rev. ST Accel. Beams 11, 114001 (2008).
  * Danilov and Nagaitsev (2010) V. Danilov and S. Nagaitsev, Phys. Rev. ST Accel. Beams 13, 084002 (2010).
  * Danilov and Nagaitsev (2014) V. Danilov and S. Nagaitsev, Phys. Rev. ST Accel. Beams 17, 124402 (2014).
  * Lichtenberg and Lieberman (1992) A. Lichtenberg and M. Lieberman,  _Regular and Chaotic Dynamics_ (Springer-Verlag New York, 1992).
  * Dragt (2013) A. Dragt, in  _Handbook of Accelerator Physics and Engineering_ , edited by A. Chao, K. Mess, M. Tigner, and F. Zimmermann (World Scientific, 2013) pp. 99–105.
  * Dilão and Alves-Pires (1996) R. Dilão and R. Alves-Pires,  _Nonlinear Dynamics in Particle Accelerators_ (World Scientific Publishing, 1996).
  * Arnold and Avez (1968) V. Arnold and A. Avez,  _Ergodic Problems of Classical Mechanics. (The Mathematical Physics Monograph Series)_ (W. A. Benjamin, Inc., New York/Amsterdam, 1968).
  * Veselov (1991) A. P. Veselov, Russian Mathematical Surveys 46, 1 (1991).
  * Meiss (1992) J. D. Meiss, Rev. Mod. Phys. 64, 795 (1992).
  * Danilov (deceased) V. Danilov (deceased), proposed this theorem in private communications (2013-02-03).
  * Iatrou and Roberts (2002) A. Iatrou and J. A. G. Roberts, Nonlinearity 15, 459 (2002).
  * Hénon (1969) M. Hénon, Quart. Appl. Math. 27, 291 (1969).
  * Dullin and Meiss (2000) H. Dullin and J. Meiss, Physica D Nonlinear Phenomena 143, 262 (2000).
  * Shiltsev  _et al._ (1999) V. Shiltsev, V. Danilov, D. Finley, and A. Sery, Phys. Rev. ST Accel. Beams 2, 071001 (1999).
  * Shiltsev  _et al._ (2008) V. Shiltsev, K. Bishofberger, V. Kamerdzhiev, S. Kozub, M. Kufer, G. Kuznetsov, A. Martinez, M. Olson, H. Pfeffer, G. Saewert, V. Scarpine, A. Seryi, N. Solyak, V. Sytnik, M. Tiunov, L. Tkachenko, D. Wildman, D. Wolff, and X.-L. Zhang, Phys. Rev. ST Accel. Beams 11, 103501 (2008).
  * Lobach  _et al._ (2018) I. Lobach, S. Nagaitsev, E. Stern, and T. Zolkin, in  _Proc. 9th International Particle Accelerator Conference (IPAC’18), Vancouver, BC, Canada, April 29-May 4, 2018_ , International Particle Accelerator Conference No. 9 (JACoW Publishing, Geneva, Switzerland, 2018) pp. 3143–3145.
  * Antipov  _et al._ (2017) S. Antipov, D. Broemmelsiek, D. Bruhwiler, D. Edstrom, E. Harms, V. Lebedev, J. Leibfritz, S. Nagaitsev, C. Park, H. Piekarz, P. Piot, E. Prebys, A. Romanov, J. Ruan, T. Sen, G. Stancari, C. Thangaraj, R. Thurman-Keup, A. Valishev, and V. Shiltsev, Journal of Instrumentation 12, T03002 (2017).
  * Brown (1983) M. Brown, The American Mathematical Monthly 90, 569 (1983).
  * Brown (1993) M. Brown, in  _Continuum Theory and Dynamical Systems (Lecture notes in pure and applied mathematics)_ , Vol. 149, edited by T. West (CRC Press, 1993) pp. 83–87.