# Dynamical structures associated with high-order and secondary resonances in
the spin-orbit problem

Hanlun Lei  School of Astronomy and Space Science, Nanjing University, Nanjing
210023, China  Key Laboratory of Modern Astronomy and Astrophysics in Ministry
of Education, Nanjing University, Nanjing 210023, China

###### Abstract

In our Solar system, spin-orbit resonances are common under Sun–planet,
planet–satellite and binary asteroid configurations. In this work, high-order
and secondary spin-orbit resonances are investigated by taking numerical and
analytical approaches. Poincaré sections as well as two types of dynamical
maps are produced, showing that there are complicated structures in the phase
space. To understand numerical structures, we adopt the theory of perturbative
treatments to formulate resonant Hamiltonian for describing spin-orbit
resonances. Results show that there is an excellent agreement between
analytical and numerical structures. It is concluded that the main V-shape
structure arising in the parameter space $(\dot{\theta},\alpha)$ is sculpted
by the synchronous primary resonance, those minute structures inside the
V-shape region are dominated by secondary resonances and those structures
outside the V-shape region are governed by high-order resonances. At last, the
analytical approach is applied to binary asteroid systems (65803) Didymos and
(4383) Suruga to reveal their phase-space structures.

celestial mechanics – minor planets, asteroids: general – planetary systems

††journal: AJ

##  1 Introduction

Spin-orbit resonance is a common phenomenon in the Universe and it takes place
if the rotational period of one body and its orbital period around the central
object are commensurable (Goldreich & Peale, 1966; Peale, 1977). For example,
the planet Mercury is located inside the 3:2 spin-orbit resonance (Lemaitre et
al., 2006; Peale & Gold, 1965), and our Moon is locked inside a synchronous
(1:1) spin-orbit resonance, resulting in the fact that its one side is
permanently facing our Earth (Peale, 1969). Additionally, chaotic rotations
are observed for minor objects in our Solar system, such as Saturn’s satellite
Hyperion (Wisdom et al., 1984), Nyx and Hydra (Correia et al., 2015; Showalter
& Hamilton, 2015). In recent years, it is found that spin-orbit resonance
commonly exists in binary asteroid systems (Pravec et al., 2016; Naidu &
Margot, 2015; Scheeres et al., 2006; Pravec et al., 2019).

The conventional model of spin-orbit resonance is formulated under the
assumption that the two bodies move around their barycenter in fixed Keplerian
orbits (Goldreich & Peale, 1966; Celletti, 1990a, b). Under the classical
model, spin-orbit resonance has been extensively explored. For Saturn’s
satellite Hyperion, Wisdom et al. (1984) studied the mechanism for the onset
of chaos through the overlap of first-order resonances (Chirikov, 1979) and
they predicted that there is a large-scale chaotic behavior if the asphericity
parameter $\alpha$ satisfies

$\alpha>\alpha_{c}=\frac{1}{2+\sqrt{14e}}.$

Under a nearly-integrable Hamiltonian model, Celletti & Chierchia (2000)
numerically explored the stability of periodic orbits related to spin-orbit
resonances by checking the eigenvalues of the monodromy matrix. Flynn & Saha
(2005) formulated second-order Hamiltonian model for describing spin-orbit
primary resonances (not secondary resonances) by taking Lie-series
transformation theory. Jafari Nadoushan & Assadian (2015, 2016) provided
explicit expressions for primary spin-orbit resonance width in terms of
asphericity and eccentricity and updated the Chirikov criterion for the onset
of chaos. In addition, they thoroughly explored the effects of asphericity and
eccentricity upon the dynamical structures by taking advantage of Poincaré
sections. Later, Nadoushan & Assadian (2016) provided a geography of varieties
of rotational resonances under a fourth-order Hamiltonian model.

If the orbit and rotation of a binary are coupled, it becomes the spin-orbit
coupling problem. In this respect, Naidu & Margot (2015) simulated the coupled
spin and orbital motions for binary asteroid systems using Poincaré surfaces
of section. Hou et al. (2017) provided explicit expressions to compute the
mutual potential, force and torque between two rigid bodies with arbitrary
shapes and mass distributions. Under the ellipsoid-ellipsoid model, Hou & Xin
(2017) presented modified Hamiltonian models for describing spin-orbit, spin-
spin and spin-orbit-spin resonances, showing that the spin-orbit resonance
center may change with the mass ratio and mutual distance and the center of
spin-spin and spin-orbit-spin resonances may change with the rotation states
of the binary asteroid. Recently, Jafari-Nadoushan (2023) offered a criterion
for dynamical closeness of asteroid pair based on the difference between
rotational and orbital angular momenta, and the author explored the dependence
of phase-space structures upon the semimajor axis, eccentricity, mass ratio of
the secondary to the primary and equatorial elongation of the secondary.

Almost all the works mentioned above concentrate on primary spin-orbit
resonances. However, if the frequency of small-amplitude oscillations about
the primary resonance and the orbital frequency are commensurable, secondary
spin-orbit resonances may happen. The existence of secondary resonances make
the phase-space structures be more complicated. In this respect, Wisdom (2004)
developed a perturbative model to study the secondary 3:1 spin-orbit resonance
of Enceladus and discussed the rate of tidal dissipation for Enceladus inside
this secondary resonance. Gkolias et al. (2016, 2019) took advantage of Lie-
series transformation theory to formulate analytical approximations for those
low-order secondary resonances (1:1, 2:1 and 3:1) bifurcating from the
synchronous primary resonance. For the bifurcation diagram of secondary 1:1
resonance, the analytical (normal form) solution can agree with numerical
result in the range of $\alpha\in[0,1.2]$, but starts to diverge when $\alpha$
is greater than 1.2 (Gkolias et al., 2019).

It is not difficult to see that the Hamiltonian model based on Lie-series
transformation adopted by Flynn & Saha (2005) can describe high-order primary
resonances outside the primary island but cannot describe secondary resonances
inside the primary island, while the Hamiltonian model based on Lie-series
method adopted by Gkolias et al. (2016, 2019) can be used to study secondary
resonances but cannot be used to study those high-order spin-orbit resonances.
Thus, it is required a unified resonant Hamiltonian model for describing high-
order and secondary resonances in the spin-orbit problem to explain phase-
space structures.

In this work, two types of dynamical maps (the fast Lyapunov indicator and the
normalized second-derivative-based index) are produced in the
$(\dot{\theta},\alpha)$ space. There are distinguishable structures, governed
by synchronous primary resonance, high-order resonances outside the primary
resonance and secondary resonances inside the primary resonance. To understand
these numerical structures, the theory of perturbative treatments (Henrard &
Lemaitre, 1986; Henrard, 1990) is taken to formulate a unified resonant
Hamiltonian model, which can be used to study both the high-order and
secondary resonances analytically. Results show that there is an excellent
correspondence between analytical structures arising in phase portraits and
numerical structures arising in dynamical maps.

The remaining part of this work is organized as follows. In Section 2, the
Hamiltonian function of the spin-orbit problem is briefly introduced. In
Section 3, two types of dynamical maps are produced and Section 4 takes
advantage of the theory of perturbative treatments to formulate the resonant
Hamiltonian model. Applications to binary asteroid systems are presented in
Section 5 and conclusions are summarized in Section 6.

##  2 Hamiltonian model

The spin-orbit problem states that a triaxial satellite moves around the
central planet on a fixed Keplerian orbit with semimajor axis $a$ and
eccentricity $e$ (Murray & Dermott, 1999). To simplify the problem, we
consider a planar configuration where the satellite’s spin-axis (along its
shortest physical axis) is assumed to be normal to its orbital plane. The
principal moments of inertia of the triaxial satellite are denoted by $A$, $B$
and $C$, satisfying $A\leq B\leq C$. For convenience, we take the orbital
semimajor axis $a$ and the mass of the central planet as the units of length
and mass, respectively, and take the time unit to make the mean motion $n$ be
unitary. In normalized units, the orbital period of satellite is equal to
$2\pi$ and the gravitational constant ${\cal G}$ becomes unitary.

The normalized Hamiltonian, governing the evolution of satellite’s spin axis,
can be written as (Goldreich & Peale, 1966; Flynn & Saha, 2005; Gkolias et
al., 2016, 2019; Celletti, 1990a, b)

${\cal
H}=\frac{1}{2}{\dot{\theta}^{2}}-\frac{{{\alpha^{2}}}}{{4{r^{3}}}}\cos\left[{2\left({\theta-f}\right)}\right],$
(1)

where $\theta$ is the angle between the satellite’s longest axis and a
reference line, $f$ is the true anomaly and $r$ is the radial distance from
the planet and the asphericity parameter $\alpha$ is characterized in terms of
the moments of inertia (or equatorial elongation of satellite) as follows
(Murray & Dermott, 1999):

$\alpha=\sqrt{\frac{{3\left({B-A}\right)}}{C}}=\sqrt{\frac{{3\left({a_{s}^{2}/b_{s}^{2}-1}\right)}}{{a_{s}^{2}/b_{s}^{2}+1}}},$

where $a_{s}$ and $b_{s}$ are the semi-axis of triaxial satellite ($\alpha$ is
independent upon the shortest semi-major axis $c_{s}$). As $r(t)$ and $f(t)$
are functions of time $t$, the Hamiltonian (1) determines a non-autonomous
dynamical model. To make it be autonomous, it is standard to augment the phase
space by introducing an action variable $T$ which is conjugated to time $t$.
The resulting new set of canonical variables are denoted by (Flynn & Saha,
2005)

$\displaystyle q_{1}=\theta,\quad p_{1}=\dot{\theta},$ $\displaystyle
q_{2}=t,\quad p_{2}=T$

and the augmented Hamiltonian becomes

${\cal
H}=\frac{1}{2}{\dot{q}_{1}^{2}}+p_{2}-\frac{{{\alpha^{2}}}}{{4{r^{3}}}}\cos\left[{2\left({q_{1}-f}\right)}\right],$
(2)

which determines a 2-DOF dynamical model. In practice, it is usual to take a
certain value of $p_{2}$ to make the Hamiltonian ${\cal
H}(q_{1},p_{1},q_{2},p_{2})$ be equal to zero.

According to elliptic expansions, the Hamiltonian (2) can be expanded in
Fourier series (Celletti, 1990a, b). As a consequence, the spin-orbit
Hamiltonian truncated at the fourth order in $e$ takes the form:

${\cal
H}=\frac{{p_{1}^{2}}}{2}+{p_{2}}-\frac{{{\alpha^{2}}}}{4}\sum\limits_{n=-6}^{2}{{{\cal
C}_{n}(e)}\cos\left({2{q_{1}}+n{q_{2}}}\right)},$ (3)

where the coefficients are related to eccentricity in the following form:

$\displaystyle{{\cal C}_{-6}}$
$\displaystyle=\frac{{533}}{{16}}{e^{4}},\;{{\cal
C}_{-5}}=\frac{{845}}{{48}}{e^{3}},\;{{\cal
C}_{-4}}=\frac{{17}}{2}{e^{2}}-\frac{{115}}{6}{e^{4}},$ $\displaystyle{{\cal
C}_{-3}}$ $\displaystyle=\frac{7}{2}e-\frac{{123}}{{16}}{e^{3}},\;{{\cal
C}_{-2}}=1-\frac{5}{2}{e^{2}}+\frac{{13}}{{16}}{e^{4}},$ $\displaystyle{{\cal
C}_{-1}}$ $\displaystyle=-\frac{1}{2}e+\frac{1}{{16}}{e^{3}},\;{{\cal
C}_{0}}=0,\;{{\cal C}_{1}}=\frac{1}{{48}}{e^{3}},\;{{\cal
C}_{2}}=\frac{1}{{24}}{e^{4}}.$

The equations of motion can be derived from Hamiltonian canonical relations
(Morbidelli, 2002),

$\displaystyle{{\dot{q}}_{1}}$ $\displaystyle=\frac{{\partial{\cal
H}}}{{\partial{p_{1}}}},\quad{{\dot{p}}_{1}}=-\frac{{\partial{\cal
H}}}{{\partial{q_{1}}}},$ (4) $\displaystyle{{\dot{q}}_{2}}$
$\displaystyle=\frac{{\partial{\cal
H}}}{{\partial{p_{2}}}},\quad{{\dot{p}}_{2}}=-\frac{{\partial{\cal
H}}}{{\partial{q_{2}}}}.$

In summary, the spin-orbit problem is ruled by two parameters: the orbital
eccentricity $e$ and the asphericity parameter $\alpha$. In Sect. 3, the
equations of motion are numerically integrated in order to produce dynamical
maps as well as Poincaré sections.

##  3 Dynamical maps

In this section, we numerically explore phase-space structures of the spin-
orbit problem by means of dynamical maps, which provides a global view about
the dynamics. In particular, two types of indicators are adopted. In producing
two types of dynamical maps, the asphericity parameter $\alpha$ and the spin
frequency $\dot{\theta}$ are distributed in the interval
$[0.2,1.5]$$\times$$[-0.5,2.5]$ with steps of $\delta\alpha=0.01$ and
$\delta\dot{\theta}=0.01$ and numerical integration is performed over 500
orbital periods.

The first one is the fast Lyapunov indicator (FLI), introduced by Froeschlé et
al. (1997), which provides a quick method for distinguishing between chaotic
and regular regions in the phase space. There are several variants for the
definition of FLI (Barrio, 2005, 2006; Skokos, 2009). In this work, we adopt
the definition of FLI, given by (Guzzo et al., 2002)

${\rm
FLI}=\mathop{\sup}\limits_{t\in\left[{{t_{0}},{t_{f}}}\right]}{\log_{10}}\frac{{\left\|{\delta{\bm{X}}\left(t\right)}\right\|}}{{\left\|{\delta{\bm{X}}\left({{t_{0}}}\right)}\right\|}}$
(5)

where $\bm{X}=(\theta,\dot{\theta})^{\rm T}$ is the state vector,
${\delta{\bm{X}}\left({{t_{0}}}\right)}$ is the initial tangent vector,
${\delta{\bm{X}}\left({{t}}\right)}$ is the deviation at the moment $t$, and
$t_{f}$ is the maximum integration time. In general, the magnitude of FLI
increases linearly with time for regular motions and it grows exponentially
with time for chaotic motions (Guzzo et al., 2002; Skokos, 2009). As a result,
the maps of FLI can clearly show regular and chaotic regions.

For the spin-orbit problem, we take the maximum integration time as 100 times
the orbital period and fix the initial angle at $\theta=0$. In practical
computation, we take the initial deviation relative to the reference orbit as
${\delta{\bm{X}}\left({{t_{0}}}\right)}=\sqrt{2}/2\times 10^{-4}$ which
include a small component along the tangent vector (our experient shows that
the choice of ${\delta{\bm{X}}\left({{t_{0}}}\right)}$ has no influence upon
the maps of FLI). Fig. 1 shows the maps of FLI in the $(\dot{\theta},\alpha)$
space for two configurations with orbital eccentricities at $e=0.01$ and
$e=0.02$. The blue color, standing for low magnitude of FLI, indicates regular
dynamics and the yellow color, corresponding to high magnitude of FLI,
indicates chaotic dynamics. From Fig. 1, we can see that (a) the structures
corresponding to $e=0.01$ and $e=0.02$ are qualitatively similar and the
latter one has wider chaotic layers than the former one, (b) the chaotic
strips in both panels presents V-shape structure in the
$(\dot{\theta},\alpha)$ space, (c) inside the V-shape region some sub-
structures caused by secondary resonances can be visibly observed and (d)
outside the V-shape region additional sub-structures caused by high-order
spin-orbit resonances are visible. As shown in the following sections, the
V-shape strips arising in the maps of FLI are due to the dynamical separatrix
of the synchronous primary resonance. Around the separatrix, chaotic layers
form the main V-shape structures.

Figure 1: Dynamical maps of the fast Lyapunov indicator (FLI) shown in the
$(\dot{\theta},\alpha)$ space for the configurations with orbital
eccentricities at $e=0.01$ (left panel) and $e=0.02$ (right panel). The index
shown in the color bar corresponds to the base 10 logarithm of FLI. In both
panels, the main V-shape structures are visible. The red dashed lines stand
for the locations of $\alpha=0.35$ and $\alpha=0.65$ and their associated
Poincaré sections are shown in Fig. 3.

The second indicator we adopt is the second-derivative-based index
$\left\|{\Delta D}\right\|$ where $D$ measures the diameter of an orbit,
developed by Daquin & Charalambous (2023). The second-derivative-based index
is a sensitive and robust indicator, which can be used to detect separatrices,
resonant webs and chaotic seas in the phase space. In particular, such an
indicator has a strong ability to distinguish minute structures in the phase
space.

For the spin-orbit problem, the action-angles variables are denoted by
$(\theta,\dot{\theta})$ and the diameter metric of an orbit is defined by

$D\left({{\theta_{0}},{\dot{\theta}_{0}}}\right)=\mathop{\max}\limits_{{t_{0}}\leq
t\leq{t_{f}}}{\dot{\theta}}\left({{\theta_{0}},{\dot{\theta}_{0}},t}\right)-\mathop{\min}\limits_{{t_{0}}\leq
t\leq{t_{f}}}{\dot{\theta}}\left({{\theta_{0}},{\dot{\theta}_{0}},t}\right)$
(6)

where $(\theta_{0},\dot{\theta}_{0})$ stand for the initial state and $t_{0}$
and $t_{f}$ are the initial and final moments of time. Based on the distance
metric $D$, a normalized second-order-derivative-based quantity is introduced
by

$\left\|{\Delta
D}\right\|=\frac{1}{D\left({{\theta_{0}},{\dot{\theta}_{0}}}\right)}\left(\left|\frac{\partial^{2}D\left({{\theta_{0}},{\dot{\theta}_{0}}}\right)}{\partial{\dot{\theta}_{0}^{2}}}\right|+\left|\frac{\partial^{2}D\left({{\theta_{0}},{\dot{\theta}_{0}}}\right)}{\partial{\theta_{0}^{2}}}\right|\right).$
(7)

It should be mentioned that the introduction of the index $\left\|{\Delta
D}\right\|$ follows from recent developments on Lagrangian descriptors and
arc-length of orbits (Daquin et al., 2022). Regarding the computation method
(central difference) about the second derivatives of the diameter, please
refer to Daquin & Charalambous (2023) for more details. Contrarily to the FLI
method, the method based on $\left\|{\Delta D}\right\|$ is free of tangent
dynamics.

The $\left\|{\Delta D}\right\|$ metric allows to sharply detect the separatrix
crossing and thus the maps of $\left\|{\Delta D}\right\|$ are able to provide
more details on minute structures in the phase space. The difference of
$\left\|{\Delta D}\right\|$ defined in this work from the conventional
definition given in Daquin & Charalambous (2023) lies in the normalization.
Our experiments show that the normalized indicator can detect those structures
with small strength.

For the spin-orbit problem, we take the maximum integration time as 100 times
the orbital period and fix the initial angle at $\theta_{0}=0$ (this setting
is the same as that made in FLI maps). Fig. 2 presents the dynamical
structures of the normalized second-derivative-based index $||\Delta D||$ in
the $(\dot{\theta},\alpha)$ space for eccentricities at $e=0.01$ and $e=0.02$.
Besides the V-shape structures, some minute structures can be visibly
detected. In particular, the structures on the left side of the V-shape
structure can be visibly observed in the maps of $||\Delta D||$, while they
are not visible in the maps of FLI. In addition, more sub-structures can be
found inside the V-shape region.

Figure 2: Dynamical maps of the second-derivative-based index shown in the
$(\dot{\theta},\alpha)$ space for the configurations with orbital
eccentricities at $e=0.01$ (left panel) and $e=0.02$ (right panel). The index
shown in the color bar corresponds to the base 10 logarithm of the normalized
second-derivative-based index $||\Delta D||$ (to make the structures be
distinguishable, the index $||\Delta D||$ is limited in the interval
[10,10000]). In both panels, the main V-shape structures are visible. Compared
to the maps of FLI shown in Fig. 1, minute structures can be visibly observed,
especially in the region on the left-hand side of the V-shape structure.
Poincaré sections corresponding to the cases of $\alpha=0.35$ and
$\alpha=0.65$ are presented in Fig. 3.

Figure 3: Nominal location of the centers of high-order and secondary
resonances determined under the Kernel Hamiltonian model together with the
Poincaré surfaces of section for spin-orbit problems with asphericity
parameters at $\alpha=0.35$ (left panel) and $\alpha=0.65$ (right panel).
Nominal locations of high-order and secondary resonances are determined by
equation (14). The eccentricity is taken as $e=0.01$. The red dashed lines
stand for the dynamical separatrix of the synchronous primary resonance,
specified by equation (10). The primary resonance happens inside the region
bounded by the dynamical separatrices.

To understand the structures arising in the dynamical maps shown in Figs. 1
and 2, we take two examples with $\alpha=0.35$ and $\alpha=0.65$, which are
explicitly marked in the left panels. To produce Poincaré sections, we record
the states every time when the variable $q_{2}$ (or $t$) satisfies

${\rm mod}\left({q_{2},2\pi}\right)={\rm mod}\left({t,2\pi}\right)=0.$

In simulation, the state variable $p_{2}$ is chosen to ensure that $\cal H$ is
equal to zero and the numerical integration is taken over 500 orbital periods.
Poincaré sections of these two representative systems are produced in the case
of $e=0.01$, as shown in Fig. 3. From the Poincaré sections, we can see that
(a) in the case of $\alpha=0.35$, the secondary 3:1 and 4:1 spin-orbit
resonances can be observed and (b) in the case of $\alpha=0.65$ only the
secondary 2:1 resonance can be found inside the synchronous primary resonance.
These secondary resonance are responsible for the structures arising inside
the V-shape region shown in Figs. 1 and 2. In both panels of Fig. 3, we can
further observe (a) chaotic layers around the separatrices of the primary
resonance, corresponding to the V-shape main structures, and (b) high-order
resonances outside the primary resonance, corresponding to the sub-structures
outside the V-shape region.

From the dynamical maps shown in Figs. 1 and 2, we may ask: how do high-order
and/or secondary spin-orbit resonances sculpt the dynamical structures arising
in the phase space? In the coming section, we aim at addressing this problem
by taking advantage of the theory of perturbative treatments developed by
Henrard & Lemaitre (1986) and Henrard (1990).

##  4 Perturbative treatments

To describe the synchronous primary resonance, we introduce the following
canonical transformation:

$\displaystyle{\sigma_{1}}$
$\displaystyle={q_{1}}-{q_{2}},\quad{\Sigma_{1}}={p_{1}},$
$\displaystyle{\sigma_{2}}$
$\displaystyle={q_{2}},\quad{\Sigma_{2}}={p_{2}}+{p_{1}},$

where the argument $\sigma_{1}$ is the resonant angle of the synchronous
primary resonance. Under such a new set of variables, the Hamiltonian can be
re-organized as follows:

$\displaystyle{\cal H}=$
$\displaystyle\frac{{\Sigma_{1}^{2}}}{2}+{{\Sigma_{2}}-{\Sigma_{1}}}-\frac{{{\alpha^{2}}}}{4}{{\cal
C}_{-2}}\left(e\right)\cos\left({2{\sigma_{1}}}\right)$ (8)
$\displaystyle-\frac{{{\alpha^{2}}}}{4}\sum\limits_{\scriptstyle-6\leq n\leq
2\hfill\atop\scriptstyle n\neq-2\hfill}{{{\cal
C}_{n}}\left(e\right)\cos\left[{2{\sigma_{1}}+\left({n+2}\right){\sigma_{2}}}\right]}.$

From the viewpoint of perturbative treatment (Henrard & Lemaitre, 1986;
Henrard, 1990), the Hamiltonian is composed of the kernel Hamiltonian
(unperturbed part) and the part of perturbation, denoted by

${\cal H}={\cal H}_{0}\left(\sigma_{1},\Sigma_{1},\Sigma_{2}\right)+{\cal
H}_{1}\left(\sigma_{1},\Sigma_{1},\sigma_{2},\Sigma_{2}\right),$ (9)

where the unperturbed part is

${{\cal
H}_{0}}=\frac{{\Sigma_{1}^{2}}}{2}+{{\Sigma_{2}}-{\Sigma_{1}}}-\frac{{{\alpha^{2}}}}{4}\cos\left({2{\sigma_{1}}}\right)$

and the perturbation part is

$\displaystyle{\cal H}_{1}=$
$\displaystyle-\frac{{{\alpha^{2}}}}{4}\left[{{\cal
C}_{-2}}\left(e\right)-1\right]\cos\left({2{\sigma_{1}}}\right)$
$\displaystyle-\frac{{{\alpha^{2}}}}{4}\sum\limits_{\scriptstyle-6\leq n\leq
2\hfill\atop\scriptstyle n\neq-2\hfill}{{{\cal
C}_{n}}\left(e\right)\cos\left[{2{\sigma_{1}}+\left({n+2}\right){\sigma_{2}}}\right]}.$

The reason why we keep those periodic terms related to $\sigma_{1}$ in the
kernel Hamiltonian is that the magnitudes of coefficients related to
$\sigma_{1}$ are at the zeroth order in $e$. In this way, the influence of
synchronous primary resonance upon high-order and/or secondary resonances are
taken into consideration.

The dynamical model specified by the kernel Hamiltonian ${\cal H}_{0}$ is of
one degree of freedom and it can be used to describe the dynamical structure
of the synchronous primary resonance. Under the perturbation of ${\cal
H}_{1}$, high-order spin-orbit resonances (outside the primary resonance) as
well as secondary resonances (inside the primary resonance) may appear in the
phase space.

In this section, we will first study the basic dynamics governed by the kernel
Hamiltonian and determine the nominal location of high-order and secondary
resonances in Sect. 4.1, then formulate resonant Hamiltonian models by taking
advantage of perturbative treatments (Henrard & Lemaitre, 1986; Henrard, 1990)
in Sect. 4.2 and results are shown in Sect. 4.3.

###  4.1 Dynamics under the Kernel Hamiltonian model

As $\sigma_{2}$ is absent from ${\cal H}_{0}$, the unperturbed Hamiltonian
model holds an integral of motion $\Sigma_{2}$. Thus, the kernel Hamiltonian
determines an integrable dynamical model. It is not difficult to get that the
stable equilibrium points are located at $(2\sigma_{1}=0,\Sigma_{1}=1)$ and
the unstable equilibrium points are located at
$(2\sigma_{1}=\pi,\Sigma_{1}=1)$. Thus, the synchronous primary resonance is
centered at $(2\sigma_{1}=0,\Sigma_{1}=1)$. The level curve of Hamiltonian
passing through the saddle point in the phase space, satisfying

$\frac{1}{2}\left(\Sigma_{1}^{2}+1\right)-{\Sigma_{1}}-\frac{{{\alpha^{2}}}}{4}\left(\cos
2{\sigma_{1}}+1\right)=0.$

As a consequence, the dynamical separatrix of the synchronous primary
resonance can be explicitly expressed by

${\Sigma_{1}}=1\pm\frac{{\sqrt{2}}}{2}\alpha\sqrt{1+\cos 2{\sigma_{1}}},$ (10)

where the symbol ‘+’ corresponds to the upper separatrix and ‘-’ corresponds
to the lower one. Usually, the distance between the upper and lower
separatrices evaluated at the resonance centre $(2\sigma_{1}=0,\Sigma_{1}=1)$
is called the resonant width. According to equation (10), the resonant width
(in normalised units) is equal to

$\Delta\Sigma_{1}=\Delta\dot{\theta}=\alpha.$

It shows that the asphericity parameter $\alpha$ stands for the resonant width
(or strength) of the synchronous primary resonance.

Dynamical separatrix plays an role in dividing the entire phase space into
regions of circulation and libration. In particular, the region bounded by the
separatrices described by equation (10) are referred to as the librating zone
and the remaining regions are called circulating zones. In addition, the level
curves of Hamiltonian in the phase space (i.e., phase portraits) stand for the
trajectories, which are exactly periodic.

Under the unperturbed Hamiltonian model, let us introduce action–angle
variables as follows (Morbidelli, 2002):

$\displaystyle\sigma_{1}^{\rm{*}}$
$\displaystyle={\sigma_{1}}-\rho\left({\sigma_{1}^{\rm{*}},\Sigma_{1}^{\rm{*}},{\Sigma_{2}^{*}}}\right),\quad\Sigma_{1}^{\rm{*}}=\frac{1}{{2\pi}}\oint{{\Sigma_{1}}{\rm
d}{\sigma_{1}}}$ (11) $\displaystyle\sigma_{2}^{*}$
$\displaystyle=\sigma_{2},\quad\Sigma_{2}^{*}=\Sigma_{2},$

where $\sigma_{1}^{\rm{*}}$ becomes a linear function of time,
$\rho\left({\sigma_{1}^{\rm{*}},\Sigma_{1}^{\rm{*}},{\Sigma_{2}^{*}}}\right)$
is a periodic function with period $T$ and $\Sigma_{1}^{\rm{*}}$ is called
Arnold action, which represents the signed area bounded by the level curves in
the phase space divided by $2\pi$ (Morbidelli, 2002). The Arnold action is
equal to zero at the centre of the synchronous resonance and it reaches the
maximum along the separatrix of the primary resonance.

The transformation given by equation (11) is canonical with the following
generating function,

$S\left({{\sigma_{1}},\Sigma_{1}^{*},\sigma_{2},\Sigma_{2}^{*}}\right)=\sigma_{2}\Sigma_{2}^{*}+\int{{\Sigma_{1}}\left({{{\cal
H}_{0}}\left({\Sigma_{1}^{*}},\Sigma_{2}^{*}\right),{\sigma_{1}}}\right){\rm
d}{\sigma_{1}}}.$

After the canonical transformation, the unperturbed Hamiltonian can be written
as a normal form,

${\cal H}_{0}\left(\sigma_{1},\Sigma_{1},\Sigma_{2}\right)\Rightarrow{{\cal
H}_{0}}\left({\Sigma_{1}^{*},\Sigma_{2}^{*}}\right)$ (12)

which determines the fundamental frequencies, given by

$\dot{\sigma}_{1}^{\rm{*}}=\frac{{\partial{{\cal
H}_{0}}\left({\Sigma_{1}^{*},\Sigma_{2}^{*}}\right)}}{{\partial\Sigma_{1}^{*}}},\quad\dot{\sigma}_{2}^{\rm{*}}=1.$
(13)

In particular, when the fundamental frequencies become commensurable, high-
order or secondary spin-orbit resonances may happen. As a result, the
condition of a certain $k_{1}$:$k_{2}$ spin-orbit resonance takes the form,

${k_{1}}\dot{\sigma}_{1}^{\rm{*}}-{k_{2}}\dot{\sigma}_{2}^{\rm{*}}=0\Rightarrow{k_{1}}\dot{\sigma}_{1}^{\rm{*}}-{k_{2}}=0.$
(14)

By solving equation (14), it is possible to determine the nominal location of
high-order and/or secondary resonances.

In Fig. 3, the nominal locations of high-order and secondary spin-orbit
resonances are shown in the phase space $(\theta,\dot{\theta})$ for the spin-
orbit problems with $\alpha=0.35$ and $\alpha=0.65$ and they are compared with
the associated Poincaré sections. In both panels of Fig. 3, the red dashed
lines stand for the dynamical separatrices of the synchronous primary
resonance, which are mathematically described by equation (10). It is observed
that around the dynamical separatrix there is a chaotic layer, which is wilder
with a higher asphericity parameter $\alpha$.

For convenience, we refer to the spin-orbit resonances inside the primary
resonance as the secondary resonances and the ones outside the primary
resonance as high-order resonances. In the case of $\alpha=0.35$, the
secondary 3:1 and 4:1 resonances can be observed and, in the case of
$\alpha=0.65$, only the secondary 2:1 resonance can be seen. Comparing the
distribution of high-order and secondary resonances to the Poincaré sections,
we can further observe that the distributions of high-order and secondary
resonances are in good agreement with the structures arising in Poincaré
sections.

Figure 4: Nominal locations of high-order and secondary spin-orbit resonances
determined under the Kernel Hamiltonian model shown in the
$(\dot{\theta},\alpha)$ space (left panel) and comparisons to dynamical maps
of the second-derivative-based index in the case of $e=0.01$ (right panel).
The map shown in the right panel is the same as that presented in the left
panel of Fig. 2. In the left panel, the shaded region stands for the parameter
space where the synchronous primary resonance happens. Those resonances inside
the primary resonance are called secondary resonances and the ones outside the
primary resonance are referred to as high-order spin-orbit resonances.

In Fig. 4, we further present the distribution of high-order and secondary
resonance in the $(\dot{\theta},\alpha)$ space and compare them to the
dynamical map of the normalized second-derivative-based index $||\Delta D||$.
The map of $||\Delta D||$ shown in the right panel is the same as the one
shown in Fig. 2. In the left panel, the red shaded region stands for the space
where the synchronous primary resonance takes place and the red dashed lines
stand for the separatrices of the primary resonance. Secondary resonances
happen inside the shaded region and high-order resonances takes place outside
the shaded region. It is observed that the nominal locations of spin-orbit
resonances are symmetric with respect to the nominal location of the primary
resonance at $\dot{\theta}=1$.

From the right panel of Fig. 4, we can observe that the distributions of high-
order and secondary resonances determined under the unperturbed Hamiltonian
model are in excellent agreement with the structures arising in the map of
$||\Delta D||$.

###  4.2 Resonant Hamiltonian model

From the previous subsection, it is known that high-order and secondary
resonances may take place in the phase space. The purpose of this section is
to formulate a Hamiltonian model for describing the resonant dynamics.

Under the transformation given by equation (11), the Hamiltonian can be re-
organized as

${\cal H}={{\cal H}_{0}}\left({\Sigma_{1}^{*},\Sigma_{2}^{*}}\right)+{{\cal
H}_{1}}\left({\sigma_{1}^{*},\sigma_{2}^{*},\Sigma_{1}^{*},\Sigma_{2}^{*}}\right),$
(15)

where the perturbation part ${\cal H}_{1}$ is much smaller than the Kernel
Hamiltonian ${\cal H}_{0}$. Therefore, it is reasonable to take ${\cal H}_{1}$
as a perturbation.

It should be mentioned that the core concept of henrard’s perturbative
treatment (Henrard & Lemaitre, 1986; Henrard, 1990) lies in transforming the
Hamiltonian form given by equation (9) to the standard form of Hamiltonian
given by equation (15). In the standard Hamiltonian, the kernel Hamiltonian is
of normal form. The standard Hamiltonian form given by equation (15) is ready
for studying resonances between $\sigma_{1}^{*}$ and $\sigma_{2}^{*}$.

In order to study the $k_{1}$:$k_{2}$ resonance between $\sigma_{1}^{*}$ and
$\sigma_{2}^{*}$, let us introduce the following transformation,

$\displaystyle{\gamma_{1}}$
$\displaystyle=\sigma_{1}^{\rm{*}}-\frac{{{k_{2}}}}{{{k_{1}}}}\sigma_{2}^{\rm{*}},\quad{\Gamma_{1}}=\Sigma_{1}^{*},$
(16) $\displaystyle{\gamma_{2}}$
$\displaystyle=\sigma_{2}^{\rm{*}},\quad{\Gamma_{2}}=\Sigma_{2}^{*}+\frac{{{k_{2}}}}{{{k_{1}}}}\Sigma_{1}^{*},$

which is canonical with the generating function

$S=\sigma_{1}^{\rm{*}}{\Gamma_{1}}+\sigma_{2}^{\rm{*}}\left({{\Gamma_{2}}-\frac{{{k_{2}}}}{{{k_{1}}}}{\Gamma_{1}}}\right).$

The resonant argument for the $k_{1}$:$k_{2}$ resonance is
${\gamma_{1}}=\sigma_{1}^{\rm{*}}-\frac{{{k_{2}}}}{{{k_{1}}}}\sigma_{2}^{\rm{*}}$.

Under the new set of canonical variables given by equation (16), the
Hamiltonian can be further expressed as

${\cal H}={{\cal H}_{0}}\left({{\Gamma_{1}},{\Gamma_{2}}}\right)+{{\cal
H}_{1}}\left({{\gamma_{1}},{\gamma_{2}},{\Gamma_{1}},{\Gamma_{2}}}\right).$
(17)

In particular, when the triaxial satellite is located inside a certain
$k_{1}$:$k_{2}$ spin-orbit resonance, the corresponding resonant argument
$\gamma_{1}=\sigma_{1}^{\rm{*}}-\frac{{{k_{2}}}}{{{k_{1}}}}\sigma_{2}^{\rm{*}}$
becomes a long-periodic angle and the other one $\gamma_{2}$ is a short-
periodic angle. It means that equation (17) determines a separable Hamiltonian
system (Henrard, 1990). As a consequence, in the long-term evolution, those
short-periodic effects can be filtered out by means of averaging theory
(corresponding to the lowest-order perturbation theory). To this end, we
perform average for the Hamiltonian over $k_{1}$ times the period of
$\gamma_{2}$, leading to the resonant Hamiltonian,

$\displaystyle{{\cal H}^{*}}$
$\displaystyle=\frac{1}{{2{k_{1}}\pi}}\int\limits_{0}^{2{k_{1}}\pi}{{\cal
H}\left({{\gamma_{1}},{\gamma_{2}},{\Gamma_{1}},{\Gamma_{2}}}\right){\rm
d}{\gamma_{2}}}$ (18) $\displaystyle={{\cal
H}_{0}}\left({{\Gamma_{1}},{\Gamma_{2}}}\right)+\frac{1}{{2{k_{1}}\pi}}\int\limits_{0}^{2{k_{1}}\pi}{{{\cal
H}_{1}}\left({{\gamma_{1}},{\gamma_{2}},{\Gamma_{1}},{\Gamma_{2}}}\right){\rm
d}{\gamma_{2}}}.$

Under the resonant Hamiltonian model determined by ${\cal H}^{*}$, the angle
$\gamma_{2}$ becomes a cyclic variable, so that its conjugate momentum

${\Gamma_{2}}=\Sigma_{2}^{*}+\frac{{{k_{2}}}}{{{k_{1}}}}\Sigma_{1}^{*}$

is a motion integral. As a result, the resonant Hamiltonian model determined
by equation (18) is integrable.

Given the motion integral, the global dynamical structures can be explored by
plotting level curves of resonant Hamiltonian ${\cal H}^{*}$ in the phase
space (i.e., phase portraits). In order to compare phase portraits with
Poincaré sections, we need to project the phase space
$(\sigma_{1}^{*},\Sigma_{1}^{*})$ to the space $(\theta,\dot{\theta})$ by
considering the condition of Poincaré sections: ${\rm
mod}\left(q_{2},2\pi\right)=0$.

Figure 5: Dynamical separatrices of high-order and secondary spin-orbit
resonances in the phase space (left-column panels) as well as their
comparisons to numerical structures arising in Poincaré sections (right-column
panels) for spin-orbit problems with different asphericity parameters. The red
dashed lines stand for dynamical separatrices of the synchronous primary
resonance, which are expressed by equation (10).

###  4.3 Results

Analytical and numerical results for high-order and secondary spin-orbit
resonances in the phase space $(\theta,\dot{\theta})$ are presented in Fig. 5.
Regarding analytical result, we provide the dynamical separatrix, which stands
for the isolines of resonant Hamiltonian. Dynamical separatrix divides the
phase space into regions of circulation and libration. In essence, dynamical
separatrix corresponds to invariant manifolds stemming from unstable periodic
orbits in the original dynamical model. In the left-column panels, dynamical
separatrices of the synchronous primary resonance, high-order resonances and
secondary resonances are shown and, in the right-column panels, they
structures are compared to numerical structures arising in Poincaré sections.
In practical simulations, the spin-orbit problems with asphericity parameters
at $\alpha=0.35$, $\alpha=0.65$, $\alpha=1.1$ and $\alpha=1.4$ are considered
and the eccentricity is taken as $e=0.01$. In all panels, the dynamical
separatrices of the synchronous primary resonance are marked in red dashed
lines.

From Fig. 5, it is observed that the analytical structures (i.e., dynamical
separatrices) of high-order and/or secondary spin-orbit resonances are in good
agreement with numerical structures arising in the Poincaré sections. In
particular, the numerical structures are dominated by the secondary 3:1 and
4:1 spin-orbit resonances and high-order resonances including 2:1, 4:1 and
2:$-1$ in the case of $\alpha=0.35$, by the secondary 2:1 resonance and high-
order resonances including 2:1 and 2:$-1$ in the case of $\alpha=0.65$ and by
the secondary 1:1 resonance and high-order 1:1 resonance in the cases of
$\alpha=1.1$ and $\alpha=1.4$.

In the spin-orbit problem with $\alpha=0.35$, the dynamical separatrices of
the secondary 4:1 resonance (inside the primary resonance) and high-order 4:1
resonance (outside the primary resonance) are close to the separatrices of the
primary resonance, and they provide boundaries for the chaotic layer around
the separatrix of the primary resonance. When the asphericity parameter is
increased to $\alpha=0.65$, the high-order 4:1 resonance disappears and the
secondary 3:1 and 4:1 resonances are replaced by the secondary 2:1 resonance.
It is observed that the dynamical separatrices of the high-order 2:1 and
secondary 2:1 resonances provide boundaries for chaotic layer around the
separatrix of the primary resonance. When the asphericity parameter is further
increased up to $\alpha=1.1$ and $\alpha=1.4$, only high-order 1:1 and
secondary 1:1 resonances can be observed in the phase space.

Figure 6: Analytical characteristic curves of the synchronous primary
resonance and the secondary 1:1 spin-orbit resonance (left panel) and
comparisons to the associated numerical curves (right panel). The vertical
dashed line stands for the nominal location of the synchronous primary
resonance and the shaded area represents the space where the synchronous
primary resonance takes place. The agreement between analytical and numerical
curves is good over the entire considered range $\alpha\in[0,1.5]$.

According to the resonant Hamiltonian model, we can analytically determine the
location of synchronous primary and/or secondary 1:1 resonances
(characteristic curves). On the other hand, the numerical characteristic
curves can be determined by computing the associated periodic orbits in a high
accuracy. Denote the period by $T$. Periodic orbits should satisfy the
following condition:

$\left\\{\begin{array}[]{l}\theta\left({t={t_{0}}}\right)=\theta\left({t={t_{0}}+T}\right),\\\
\dot{\theta}\left({t={t_{0}}}\right)=\dot{\theta}\left({t={t_{0}}+T}\right).\end{array}\right.$

Without loss of generality, we assume the angle at $t=t_{0}$ is equal to zero,
i.e. ${\theta}_{0}=0$. Thus, there are only two variables
$({\dot{\theta}}_{0},T)$ to be determined for a periodic orbit. In practice,
we adopt the Newton-Raphson method to solve the unknown variables of periodic
orbits.

Fig. 6 shows the analytical characteristic curves in the left panel and
comparisons between analytical and numerical curves in the right panel. The
shaded region stands for the space where the synchronous primary resonance
happens. From Fig. 6, we can see that (a) the analytical characteristic curves
can agree well with the numerical ones in the entire considered range
$\alpha\in[0,1.5]$, and (b) the secondary 1:1 resonance appears when the
asphericity parameter is higher than $\alpha$$\sim$$1.06$ in the case of
$e=0.01$.

Regarding spin-orbit problem, Gkolias et al. (2016, 2019) took advantage of
Lie-series transformation theory to study the low-order secondary resonances
bifurcating from the synchronous primary resonance. Analytical characteristic
curves are produced from the 11th normal form construction as well as from
wisdom formula derived from nonlinear method of Bogoliubov and Mitropolsky
(Wisdom et al., 1984), and analytical results are compared to numerical
computations. It is concluded that the normal form method has a significant
improvement compared to the nonlinear method of Bogoliubov and Mitropolsky. In
particular, the normal form solution can agree with numerical result in the
range of $\alpha\in[0,1.2]$, but starts to diverge when $\alpha$ is greater
than 1.2 (please refer to Fig. 2 in their work for more details). However, our
analytical estimate can match the numerical solution in the entire considered
range of $\alpha$ up to 1.5 (see the right panel of Fig. 6).

Fig. 7 shows all the analytical characteristic curves of the synchronous
primary resonance, secondary resonances and high-order resonances in the left
panel and makes a comparison between analytical curves to the numerical map of
$||\Delta D||$ in the right panel. We can see that the structures arising in
the dynamical map have a good correspondence to the analytical characteristic
curves. Thus, these numerical structures can be well understood thanks to the
associated spin-orbit resonance dynamics.

In particular, Fig. 7 provides a global view for the rotational dynamics of
the secondary in the $(\dot{\theta},\alpha)$ space. It is possible for us to
predict possible dynamical outcomes for binary asteroid systems according to
the associated asphericity parameters. For example, the following predictions
can be made: (a) secondary 1:1 resonance may appear when $\alpha>1.06$, (b)
high-order 2:1 and 2:$-1$ resonances as well as secondary 2:1, 3:1, and 4:1
resonances may disappear in the high-$\alpha$ space because they are inside
the chaotic sea around the separatrix of the synchronous primary resonance,
and (c) high-order 1:1 and 1:$-1$ resonances can appear in the entire
considered space but the strength of high-order 1:$-1$ resonance is relatively
weak.

Table 1: Physical parameters for the example binary asteroid systems adopted in this work (please refer to http://www.asu.cas.cz/asteroid/binastdata.htm for detailed data of binary asteroid systems). The ratio of spin angular momentum to the orbital angular momentum for the secondary $L_{\rm spin}/L_{\rm orb}$ is estimated from the method given in Jafari-Nadoushan (2023). binary system | $D_{s}/D_{p}$ | $m_{s}/m_{p}$ |  $P_{p}$ (h) |  $P_{\rm orb}$ (h) |  $P_{s}$ (h) | $a_{s}/b_{s}$ | $a_{\rm orb}/D_{s}$ | $e_{\rm orb}$ | $L_{\rm spin}/L_{\rm orb}$ | $\alpha$
---|---|---|---|---|---|---|---|---|---|---
(65803) Didymos | 0.227 | 0.012 | 2.2593 | 11.91 | 11.91 | 1.55 | 6.706 | $<0.03$ | $2.2\times 10^{-3}$ | 1.1
(4383) Suruga | 0.188 | 0.007 | 3.4068 | 16.34 | 16.34 | 2.04 | 10.0 | $-$ | $9.9\times 10^{-4}$ | 1.36

Figure 7: Analytical characteristic curves of the synchronous primary
resonance, high-order resonances (outside the primary resonance) and secondary
resonances (inside the primary resonance) in the $(\dot{\theta},\alpha)$ space
(left panel) and comparisons with dynamical maps of the normalized second-
derivative-based index $||\Delta D||$ (right panel). In the right panel, the
locations of two typical binary-asteroid systems are marked. The asphericity
parameter is $\alpha=1.1$ for the binary-asteroid system (65803) Didymos and
it is 1.36 for (4383) Suruga. The vertical dashed line stands for the nominal
location of the synchronous primary resonance and the shaded region in the
left panel represents the space where the synchronous primary resonance
happens.

##  5 Applications to binary asteroid systems

The analytical approach is applied to two representative binary asteroid
systems: (65803) Didymos and (4383) Suruga (Pravec et al., 2016; Warner,
2013). (65803) Didymos is the target of the DART mission (Agrusa et al.,
2022). It should be mentioned that the analytical method presented in this
work is also applicable for other binary asteroid systems. The location of
these two systems are marked in the right panel of Fig. 7, which shows that
their asphericity parameters are greater than the critical value
$\alpha_{c}=1.06$ at which the secondary 1:1 resonance bifurcates from the
synchronous primary resonance. The physical parameters of the adopted binary
asteroid systems are shown in Table 1, including the ratio of the secondary
diameter to the primary diameter $D_{s}/D_{p}$, ratio of the secondary mass to
the primary mass $m_{s}/m_{p}$, spin period of the primary $P_{p}$, mutual
orbital period of system $P_{\rm orb}$, spin period of the secondary $P_{s}$,
the equatorial elongation of the secondary $a_{s}/b_{s}$, ratio of the orbital
semimajor axis to the secondary’s diameter $a_{\rm orb}/D_{s}$, mutual orbital
eccentricity $e_{\rm orb}$, ratio of spin angular momentum to the orbital
angular momentum $L_{\rm spin}/L_{\rm orb}$ and asphericity parameter
$\alpha$.

Regarding binary asteroid systems, Jafari-Nadoushan (2023) offered a criterion
for the definition of dynamical closeness of asteroid pairs: a binary asteroid
system is a close pair if the ratio of rotational angular momentum to the
orbital angular momentum is greater than 10 percent. According to their
estimate, we can get the ratios of spin angular momentum to the orbital
angular momentum are $L_{\rm spin}/L_{\rm orb}=2.2\times 10^{-3}$ for (65803)
Didymos and $L_{\rm spin}/L_{\rm orb}=9.9\times 10^{-4}$ for (4383) Suruga,
which are much smaller than 10 percent. Thus, both binary asteroid systems are
far away from close binaries, which means that the couple between orbital and
rotational evolution is very weak. In the dynamical evolution, the influence
of rotational evolution upon the orbital dynamics can be ignored. As a result,
the spin-orbit problem discussed in this work is a good approximation for such
a kind of weakly-coupled binary systems.

According to Table 1, the asphericity parameter is $\alpha=1.1$ for binary
asteroid system (65803) Didymos and it is equal to $\alpha=1.36$ for binary
asteroid system (4383) Suruga. Considering that the orbital eccentricities are
small and highly uncertain, we take two small eccentricities at $e=0.01$ and
$e=0.02$ in practical simulations.

Analytical and numerical results are presented in Fig. 8, where the top panels
are for (65803) Didymos and the bottom panels are for (4383) Suruga. The
dynamical separatrix of the synchronous primary resonance is shown in red
dashed line and the separatrix of the secondary 1:1 resonance is given in blue
line. From Fig. 8, we can observe that (a) there is an excellent agreement
between the analytical and numerical structures, and (b) chaotic layer becomes
wilder with a higher eccentricity.

According to the current physical parameters, it is highly possible that the
binary asteroid systems (65803) Didymos and (4383) Suruga are inside the
secondary 1:1 spin-orbit resonance. However, there are a large uncertainty for
their elongation $a_{s}/b_{s}$, leading to a large uncertainty for determining
asphericity parameter $\alpha$. For example, it will be outside the secondary
1:1 resonance if the uncertainty of asphericity parameter makes $\alpha$ be
smaller than $\sim$1.06. Thus, more observations are required to further
determine the asphericity parameter $\alpha$ accurately in order to predict
their real resonant states. However, Fig. 7 can help to predict other possible
resonant states when $\alpha$ is varied.

Figure 8: Analytical structures (dynamical separatrices) of the synchronous
primary resonance and the secondary 1:1 spin-orbit resonance together with the
associated numerical structures arising in Poincaré sections at different
eccentricities for binary asteroid systems (65803) Didymos (top panels) and
(4383) Suruga (bottom panels). The blue lines stand for the dynamical
separatrices of the 1:1 secondary resonance and red dashed lines stand for the
dynamical separatrices of the synchronous primary resonance. The asphericity
parameter is $\alpha=1.1$ for binary asteroid system (65803) Didymos and it is
$\alpha=1.36$ for binary asteroid system (4383) Suruga.

##  6 Conclusions

In this work, high-order and secondary resonances in the spin-orbit problem
are investigated. Dynamical maps are numerically produced in order to provide
a global picture for the dynamics. In particular, two types of dynamical
indicators are introduced. The first one is the fast Lyapunov indicator (FLI)
and the second one is the normalized second-derivative-based index $||\Delta
D||$. The maps of FLI can distinguish the regular and chaotic regions and the
maps of $||\Delta D||$ has an ability to detect minute structures in the phase
space. It is concluded that the main V-shape structures are sculpted by the
synchronous primary resonance, those minute structures inside the V-shape
region are governed by secondary resonances and the ones outside the V-shape
region are dominated by high-order spin-orbit resonances.

The main purpose of this work is to understand dynamical causes of numerical
structures. To this end, the theory of perturbative treatments is adopted. In
particular, the Hamiltonian is divided into the kernel Hamiltonian and the
perturbation part. The terms related to the synchronous primary resonance are
included in the kernel Hamiltonian, so that it is possible to consider the
influence of the synchronous primary resonance upon high-order and/or
secondary resonances. Under the kernel Hamiltonian model, a series of
canonical transformations are introduced and nominal locations of high-order
and secondary resonances are determined. By taking advantage of Henrard’s
perturbative treatments (Henrard, 1990), resonant Hamiltonian model is
formulated. Thus, it is ready to study high-order and secondary resonance
dynamics analytically. Results show that there is an excellent correspondence
between analytical structures arising in phase portraits and numerical
structures arising both in the dynamical maps and in the Poincaré sections.
Applications to two binary asteroid systems (including Didymos and Suruga)
show that these two binary systems are possible to be inside the secondary 1:1
spin-orbit resonance.

The main results are presented in Fig. 7 where the distributions of
synchronous primary resonance, high-order resonance and secondary resonances
are presented and compared with dynamical map of $||\Delta D||$. It provides a
global picture for us to predict possible resonant states when the asphericity
parameter is changed.

This work is supported by the National Natural Science Foundation of China
(Nos. 12073011, and 12233003) and the National Key R&D Program of China (No.
2019YFA0706601).

## References

  * Agrusa et al. (2022) Agrusa, H. F., Ferrari, F., Zhang, Y., Richardson, D. C., & Michel, P. 2022, The Planetary Science Journal, 3, 158
  * Barrio (2005) Barrio, R. 2005, Chaos, Solitons & Fractals, 25, 711
  * Barrio (2006) —. 2006, International Journal of Bifurcation and Chaos, 16, 2777
  * Celletti (1990a) Celletti, A. 1990a, Z. Ang. Math. Phys., 41, 174
  * Celletti (1990b) —. 1990b, Z. Ang. Math. Phys., 41, 453
  * Celletti & Chierchia (2000) Celletti, A., & Chierchia, L. 2000, CeMDA, 76, 229
  * Chirikov (1979) Chirikov, B. V. 1979, Physics reports, 52, 263
  * Correia et al. (2015) Correia, A. C., Leleu, A., Rambaux, N., & Robutel, P. 2015, Astronomy & Astrophysics, 580, L14
  * Daquin & Charalambous (2023) Daquin, J., & Charalambous, C. 2023, CeMDA, 135, 31
  * Daquin et al. (2022) Daquin, J., Pédenon-Orlanducci, R., Agaoglou, M., García-Sánchez, G., & Mancho, A. M. 2022, Physica D: Nonlinear Phenomena, 442, 133520
  * Flynn & Saha (2005) Flynn, A. E., & Saha, P. 2005, AJ, 130, 295
  * Froeschlé et al. (1997) Froeschlé, C., Lega, E., & Gonczi, R. 1997, Celestial Mechanics and Dynamical Astronomy, 67, 41
  * Gkolias et al. (2016) Gkolias, I., Celletti, A., Efthymiopoulos, C., & Pucacco, G. 2016, MNRAS, 459, 1327
  * Gkolias et al. (2019) Gkolias, I., Efthymiopoulos, C., Celletti, A., & Pucacco, G. 2019, Communications in Nonlinear Science and Numerical Simulation, 77, 181
  * Goldreich & Peale (1966) Goldreich, P., & Peale, S. 1966, AJ, 71, 425
  * Guzzo et al. (2002) Guzzo, M., Lega, E., & Froeschlé, C. 2002, Physica D: Nonlinear Phenomena, 163, 1
  * Henrard (1990) Henrard, J. 1990, CeMDA, 49, 43
  * Henrard & Lemaitre (1986) Henrard, J., & Lemaitre, A. 1986, Celest. Mech., 39, 213
  * Hou et al. (2017) Hou, X., Scheeres, D. J., & Xin, X. 2017, Celestial Mechanics and Dynamical Astronomy, 127, 369
  * Hou & Xin (2017) Hou, X., & Xin, X. 2017, AJ, 154, 257
  * Jafari-Nadoushan (2023) Jafari-Nadoushan, M. 2023, MNRAS, 520, 3514
  * Jafari Nadoushan & Assadian (2015) Jafari Nadoushan, M., & Assadian, N. 2015, Nonlinear Dynamics, 81, 2031
  * Jafari Nadoushan & Assadian (2016) —. 2016, Nonlinear Dynamics, 85, 1837
  * Lemaitre et al. (2006) Lemaitre, A., D’Hoedt, S., & Rambaux, N. 2006, Celestial Mechanics and Dynamical Astronomy, 95, 213
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics: aspects of solar system dynamics (Taylor & Francis, London and New York)
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Nadoushan & Assadian (2016) Nadoushan, M. J., & Assadian, N. 2016, Icarus, 265, 175
  * Naidu & Margot (2015) Naidu, S. P., & Margot, J.-L. 2015, The Astronomical Journal, 149, 80
  * Peale & Gold (1965) Peale, S., & Gold, T. 1965, Nature, 206, 1240
  * Peale (1969) Peale, S. J. 1969, The Astronomical Journal, 74, 483
  * Peale (1977) Peale, S. J. 1977, in Iau colloq. 28: Planetary satellites, 87
  * Pravec et al. (2016) Pravec, P., Scheirich, P., Kušnirák, P., et al. 2016, Icarus, 267, 267
  * Pravec et al. (2019) Pravec, P., Fatka, P., Vokrouhlickỳ, D., et al. 2019, Icarus, 333, 429
  * Scheeres et al. (2006) Scheeres, D. J., Fahnestock, E. G., Ostro, S. J., et al. 2006, Science, 314, 1280
  * Showalter & Hamilton (2015) Showalter, M., & Hamilton, D. 2015, Nature, 522, 45
  * Skokos (2009) Skokos, C. 2009, in Dynamics of Small Solar System Bodies and Exoplanets (Springer), 63–135
  * Warner (2013) Warner, B. D. 2013, The Minor planet bulletin, 40, 119
  * Wisdom (2004) Wisdom, J. 2004, AJ, 128, 484
  * Wisdom et al. (1984) Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137


# Metrological robustness of high photon number optical cat states

Tomás Fernández Martos  ICFO - Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Maciej Lewenstein  ICFO - Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
ICREA, Pg. Lluis Companys 23, 08010 Barcelona, Spain  Grzegorz Rajchel-
Mieldzioć<EMAIL_ADDRESS>ICFO - Institut de Ciències Fotòniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Philipp Stammer  ICFO - Institut de Ciències Fotòniques,
The Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain

###### Abstract

In the domain of quantum metrology cat states have demonstrated their utility
despite their inherent fragility with respect to photon loss. Here, we
introduce noise robust optical cat states which exhibit a metrological
robustness for phase estimation in the regime of high photon numbers. These
cat states are obtained from the intense laser driven process of high harmonic
generation (HHG), and in the ideal case of vanishing losses, show almost twice
the quantum Fisher information (QFI) compared to the even and odd cat states.
However, and more importantly, these HHG-cat states are much more robust
against noise such that the noisy HHG-cat outperforms the pure even/odd cat
states even in the presence of more than $25\%$ losses in the regime of high
photon numbers. Furthermore, in the regime of small losses, the HHG-cat
remains almost pure while the even/odd cat state counterpart already decohere
to the maximally mixed state. This demonstrates that high photon number
optical cat states can indeed be used for metrological applications even in
the presence of losses.

##  I Introduction

Quantum metrology exploits properties such as entanglement or squeezing to
achieve highly sensitive measurements of physical parameters with potentially
higher precision than can be achieved using purely classical approaches
giovannetti2004quantum ; giovannetti2011advances . A standard way of
quantifying the metrological limits for parameter estimation is the use of the
quantum Fisher information (QFI) liu2020quantum . The QFI sets an upper bound
on the achievable information gain from any measurement of a quantum state
$\hat{\rho}(\bm{\theta})$, with encoded unknown parameters $\bm{\theta}$,
under the optimal choice of observables. This is directly tied to a lower
bound on the achievable precision of the measurement, given by the quantum
Cramér-Rao bound helstrom1969quantum . The QFI gives a measure of the
sensitivity of the quantum state $\hat{\rho}(\bm{\theta})$ to a change of a
certain parameter from the set $\bm{\theta}$; by maximizing this quantity,
upper bounds on the sensitivity for specific measurements can be found.
However, for any realistic setting one needs to take into account the effect
of noise, which limits the attainable precision of measurements
toth2014quantum ; meyer2021fisher .

One of the possibilities for surpassing classical limits of metrology is to
use quantum superpositions of coherent states, known as Schrödinger cat states
gerry1997quantum ; gilchrist2004schrodinger ; ralph2003quantum , which can be
advantageous for phase estimation joo2011quantum ; lee2015quantum . Cat states
can be created in various systems, e.g., employing one-axis twisting dynamics
Kitagawa_1993 ; Wineland_1994 , which can be implemented with ultra-cold atoms
in optical lattices Plodzien2020 ; Plodzien2022 ; Plodzien2023generation ;
Hernandez2022 ; Dziurawiec2023 ; yanes2023spin . Here, we focus on optical cat
states, conventionally given by the superposition of two coherent states of
opposite phase, the most archetypal being the even/odd cat states
gerry1997quantum . The coherent states in the superposition are
macroscopically distinguishable for high values of the mean photon number,
i.e. coherent state amplitude. In particular, for phase estimation the
measurement uncertainty is typically inversely proportional to the average
photon number escher2011general . However, to savor the full potential of this
scaling high photon number optical cat states are required, while conventional
approaches are limited by their low photon number ourjoumtsev2006generating ;
ourjoumtsev2007generation . In addition, cat states are very vulnerable to the
decoherence introduced by photon loss for increasing photon number
glancy2008methods ; zhang2013quantum . Nonetheless, optical cat states are
still used in quantum metrology since they often provide an enhanced phase
sensitivity yu2018maximal .

In addition to the conventional approaches for generating optical cat states
ourjoumtsev2007generation ; glancy2008methods , there have been recent
advances in the generation of high photon number optical cat states, by using
the platform of intense laser matter interaction for quantum state engineering
of light lewenstein2021generation ; stammer2023quantum ;
bhattacharya2023strong . In particular, the processes of high harmonic
generation (HHG) lewenstein2021generation ; rivera2022strong  and above-
threshold ionization rivera2022light  have been utilized to generate non-
classical states of light, especially massive entangled states stammer2022high
; stammer2022theory ; stammer2023entanglement  and high photon number optical
cat states stammer2022high ; rivera2022strong ; rivera2021new . The resulting
HHG-cat state is given by the superposition of two coherent states with
amplitude $\alpha$, which have a relative amplitude displacement of
$\delta\alpha$. By tuning the relative displacement $\delta\alpha$, the
macroscopic distinguishability between the two coherent states can be
controlled. In theory, the mean photon number of the HHG-cat state can be made
arbitrarily large, by increasing $\alpha$, without affecting their resilience
against decoherence (which depends on $\delta\alpha$). This makes them
potentially more suitable for quantum metrology applications in realistic
scenarios including noise than conventional cat states. In practice, the
generation of an optical HHG-cat state with high average photon number of
$\expectationvalue{N}\approx 150$ has been shown in Ref. lamprou2023nonlinear
, orders of magnitude higher than those of conventional schemes. In this work,
we want to explore the usefulness of the optical cat state obtained from HHG
for metrological applications by analyzing its QFI, as well as to show its
resilience against photon loss in the high photon number regime. This
ultimately allows to make use of the features arising from the quantum
superposition in the high photon number regime.

##  II The ideal scenario

The aforementioned HHG-cat state is obtained when driving atoms by an intense
laser field $\ket{\alpha}$ and performing conditioning approaches on the
process of HHG lewenstein2021generation . These conditioning measurements,
extensively described in stammer2022high ; stammer2022theory ;
stammer2023quantum , lead to a quantum operation acting on the driving field
mode and result in the generation of a high photon number optical cat state.
It has been shown by tomographic reconstruction lewenstein2021generation ;
rivera2022strong  of the state that it has the following form

$\ket{\psi_{\mathrm{HHG}}}=\frac{1}{\sqrt{N_{\mathrm{HHG}}}}\left(\ket{\alpha+\delta\alpha}-\xi\ket{\alpha}\right),$
(1)

where $\xi=\bra{\alpha}\ket{\alpha+\delta\alpha}$ and
$N_{\mathrm{HHG}}=1-\absolutevalue{\xi}^{2}$ is the normalization. This
approach can be used for quantum state engineering of light stammer2023quantum
; stammer2023role ; stammer2023entanglement  in different target systems
rivera2022quantum ; rivera2023bloch  or when considering driving the processes
with non-classical states of light gorlach2023high ; stammer2023limitations ;
pizzi2023light .

To demonstrate the metrological robustness of the optical cat state generated
form HHG we consider the problem of phase estimation. The setup under
consideration is a beam splitter (BS) parameterized by the angle $\theta$
campos1989quantum , and described by the unitary operation
$\hat{U}_{\mathrm{BS}}(\theta)=e^{-i\theta\hat{H}}$, with the BS Hamiltonian
$\hat{H}=-i(\hat{a}_{0}^{\dagger}\hat{a}_{1}-\hat{a}_{0}\hat{a}_{1}^{\dagger})$.
In particular, we consider a homodyne configuration in which the input signal
state $\hat{\rho}_{0}$ is overlapped with a local oscillator (LO) in a
coherent state $\ket{\beta_{1}}$. The output state after the BS is then given
by
$\hat{\rho}_{\mathrm{out}}=\hat{U}_{\mathrm{BS}}(\theta)\hat{\rho}_{\mathrm{in}}\hat{U}^{\dagger}_{\mathrm{BS}}(\theta)$,
where
$\hat{\rho}_{\mathrm{in}}=\hat{\rho}_{0}\otimes\outerproduct{\beta_{1}}{\beta_{1}}$.
For the single parameter BS and arbitrary signal state $\hat{\rho}_{0}$, the
quantum Fisher information (QFI) can be written as yu2018maximal

$\begin{split}\mathcal{F}=&\sum_{i=1}^{d}4\lambda_{i}\bra{\lambda_{i}}\hat{H}^{2}\ket{\lambda_{i}}\\\
&-\sum_{i,j=1}^{d}\frac{8\lambda_{i}\lambda_{j}}{\lambda_{i}+\lambda_{j}}|\bra{\lambda_{i}}\hat{H}\ket{\lambda_{j}}|^{2},\end{split}$
(2)

where
$\hat{\rho}_{\mathrm{in}}=\sum_{i=1}^{d}\lambda_{i}\outerproduct{\lambda_{i}}{\lambda_{i}}$
represents the spectral decomposition of the density matrix for the input
state, while $d$ is the dimension of the matrix support, $\lambda_{i}$ are the
eigenvalues, and $\ket{\lambda_{i}}$ are the eigenvectors. For a pure input
state,
$\hat{\rho}_{\mathrm{in}}=\outerproduct*{\Psi_{\mathrm{in}}}{\Psi_{\mathrm{in}}}$,
this expression simplifies to

$\mathcal{F}=4\left(\bra{\Psi_{\mathrm{in}}}\hat{H}^{2}\ket{\Psi_{\mathrm{in}}}-\bra{\Psi_{\mathrm{in}}}\hat{H}\ket{\Psi_{\mathrm{in}}}^{2}\right).$
(3)

First, we consider the ideal scenario in which neither losses nor decoherence
affect the optical cat states. In particular, our aim is to compute the QFI
for the HHG-cat state, and show its advantage over the traditional even $(+)$
and odd $(-)$ cat states, which are superpositions of coherent states with
opposite phase
$\ket{\psi_{\pm}}=(N_{\pm})^{-1/2}(\ket{\alpha}\pm\ket{-\alpha})$, where
$N_{\pm}$ are the normalization factors gerry1997quantum . For all input
states, the first term in Eq. (3) reads

$\begin{split}\langle\Psi_{\mathrm{in}}|\hat{H}^{2}|\Psi_{\mathrm{in}}\rangle=&-2\Re{(\beta^{*})^{2}\langle\hat{a}_{0}^{2}\rangle}+|\beta|^{2}\\\
&+\langle\hat{N}_{0}\rangle(1+2|\beta|^{2}),\end{split}$ (4)

where $\hat{N}_{0}=\hat{a}^{\dagger}_{0}\hat{a}_{0}$ is the photon number
operator for the signal input mode. For the second term in Eq. (3), we obtain
$\langle\hat{H}\rangle_{\mathrm{HHG}}=2\Im{(\alpha+\delta\alpha)\beta^{*}}$
for the HHG-cat state and $\langle\hat{H}\rangle_{\pm}=0$ vanishes for the
even/odd cat state. In Fig. 1 we show the advantage of the HHG-cat over the
even/odd cat states for a wide range of photon numbers. We observe that the
HHG-cat state outperform the even/odd counterparts in the entire region of
high photon numbers $|\alpha|^{2}\geq 4$ considered in this work, and as the
coherent state amplitude increases, the enhancement of the QFI stabilizes at
around $\mathcal{F}_{\mathrm{HHG}}\approx 1.88\times\mathcal{F}_{\pm}$.
However, we shall now consider the realistic scenario in which the pure cat
states suffer from losses and decoherence, and analyze the metrological
robustness of the HHG-cat state over the even/odd cat states.

Figure 1: QFI for the phase sensing setup using different pure signal states
overlapped with a coherent state LO, normalized to the QFI for a pure odd cat
state ($\mathcal{F}/\mathcal{F}_{-}$). The coherent state amplitude of the
even/odd cat states is given by $\absolutevalue{\alpha}$, and for the HHG-cat
state we have $\alpha=\alpha^{\prime}+\delta\alpha$ for a fixed displacement
$\delta\alpha=-0.5$, such that $\alpha^{\prime}\in[2.5,10.5]$ (therefore the
mean photon numbers of all input states only differ by less than a single
photon). The mean photon number of the LO is set to match the mean photon
number of the input cat states, $|\beta|^{2}=|\alpha|^{2}$. For high values of
the mean photon number, the QFI for the HHG-cat state input converges towards
$\mathcal{F}_{\mathrm{HHG}}\approx 1.88\times\mathcal{F}_{\pm}$.

##  III Noise robust optical cat states

Typically, the primary source of decoherence in optical experiments comes from
photon losses, which can be modeled by introducing a fictitious BS with
transmittance $\eta$ and tracing out the reflected part of the output
glancy2008methods ; zhang2021quantifying ; rohde2007practical .
Mathematically, this can be expressed as

$\hat{\rho}_{\mathrm{out},\mathrm{t}}^{\eta}=\Tr_{\mathrm{r}}\big{[}\outerproduct*{\Psi_{\mathrm{out},\mathrm{tr}}^{\eta}}{\Psi_{\mathrm{out},\mathrm{tr}}^{\eta}}\big{]},$
(5)

where
$\ket*{\Psi_{\mathrm{out},\mathrm{tr}}^{\eta}}=\hat{U}_{\mathrm{BS}}\big{[}\\!\arccos{(\sqrt{\eta})}\big{]}\ket{\psi}_{0}\ket{0}_{1}$.
For the even/odd-cat states we have the mixture

$\begin{split}\hat{\rho}_{\pm}^{\eta}=(1-\zeta_{\pm}^{\eta})\outerproduct*{\psi^{\eta}_{\pm}}{\psi^{\eta}_{\pm}}+\zeta_{\pm}^{\eta}\outerproduct*{\psi^{\eta}_{\mp}}{\psi^{\eta}_{\mp}},\end{split}$
(6)

where
$\ket*{\psi^{\eta}_{\pm}}=(N_{\pm}^{\eta})^{-1/2}\big{(}\ket*{\alpha\sqrt{\eta}}\pm\ket*{-\alpha\sqrt{\eta}}\big{)}$
are the normalized even/odd cat states with a reduced amplitude by the factor
$\sqrt{\eta}$, and the parameter that determines the purity of
$\hat{\rho}_{\pm}^{\eta}$ is
$\zeta_{\pm}^{\eta}=(N_{\mp}^{\eta}/2N_{\pm})(1-\exp[-2|\alpha|^{2}(1-\eta)])$.
The purity of the mixture
$\gamma_{\pm}^{\eta}\equiv\Tr{(\hat{\rho}_{\pm}^{\eta})^{2}}=(\zeta_{\pm}^{\eta})^{2}+(1-\zeta_{\pm}^{\eta})^{2}$
rapidly decreases for increasing photon number even in the presence of small
losses serafini2004minimum . Therefore, for any realistic scenario the
even/odd cat states will rapidly decohere into a 50/50 mixture of even and odd
cat states with the reduced amplitude $\alpha\sqrt{\eta}$ due to losses. We
now consider the pure HHG-cat state in Eq. (1), and the resulting mixed state
in its spectral decomposition, after the noise channel, reads

$\hat{\rho}_{\mathrm{HHG}}^{\eta}=\lambda^{\eta}_{(+)}\outerproduct*{\lambda^{\eta}_{(+)}}{\lambda^{\eta}_{(+)}}+\lambda^{\eta}_{(-)}\outerproduct*{\lambda^{\eta}_{(-)}}{\lambda^{\eta}_{(-)}},$
(7)

where the eigenvalues are given by

$\begin{split}\lambda^{\eta}_{(\pm)}=&\frac{1}{2N_{\mathrm{HHG}}}\Big{[}N_{\mathrm{HHG}}\pm\big{[}N_{\mathrm{HHG}}^{2}\\\
&-4N^{\eta}_{\mathrm{HHG}}(|\xi^{\eta}|^{2}-|\xi|^{2}-|\xi^{\eta}-\xi\mu^{\eta}|^{2})\big{]}^{1/2}\Big{]},\end{split}$
(8)

with
$\xi^{\eta}=\innerproduct{\alpha\sqrt{\eta}}{(\alpha+\delta\alpha)\sqrt{\eta}}$,
$N^{\eta}_{\mathrm{HHG}}=1-|\xi^{\eta}|^{2}$, and
$\mu^{\eta}=\innerproduct*{(\alpha+\delta\alpha)\sqrt{1-\eta}}{\alpha\sqrt{1-\eta}}$,
and with the corresponding eigenvectors

$\begin{split}\ket*{\lambda^{\eta}_{(\pm)}}=&\frac{1}{\sqrt{N^{\eta}_{(\pm)}}}\big{[}C_{(\pm),1}^{\eta}\ket*{(\alpha+\delta\alpha)\sqrt{\eta}}\\\
&+C_{(\pm),2}^{\eta}\ket*{\alpha\sqrt{\eta}}\big{]},\end{split}$ (9)

where $N^{\eta}_{(\pm)}$ are normalization factors, and

$\begin{split}C_{(\pm),1}^{\eta}=&\frac{((\xi^{\eta})^{*}-\xi^{*}(\mu^{\eta})^{*})\sqrt{N^{\eta}_{\mathrm{HHG}}}}{N_{\mathrm{HHG}}\lambda^{\eta}_{(\pm)}-N^{\eta}_{\mathrm{HHG}}},\\\
C_{(\pm),2}^{\eta}=&\sqrt{N^{\eta}_{\mathrm{HHG}}}-\xi^{\eta}C_{(\pm),1}^{\eta}.\end{split}$
(10)

In this case, the purity of the mixture is given by
$\gamma_{\mathrm{HHG}}^{\eta}=(\lambda^{\eta}_{(+)})^{2}+(\lambda^{\eta}_{(-)})^{2}$,
and we show the purity of the different cat states as a function of $\eta$ in
Fig. 2 (a) for an initial photon number of $\langle\hat{N}\rangle=100$. We
clearly see the rapid decrease of the purity to the maximally mixed state with
$\gamma_{\pm}^{\eta}=1/2$ for the even/odd cat state in this photon number
regime already for very small losses serafini2004minimum . In contrast, the
HHG-cat state never becomes maximally mixed, i.e.
$\gamma_{\mathrm{HHG}}^{\eta}>1/2$ for the 2-dimensional decomposition, and
outperforms the even/odd cat states for all $\eta$. Notably, the purity of the
HHG-cat is not affected by the average photon number since the eigenvalues in
(8) only depend on the relative displacement $\delta\alpha$, while in
contrast, the even/odd cat states asymptotically tend to a constant
$\gamma_{\pm}^{\eta}=1/2$ for increasing $|\alpha|$ serafini2004minimum .
Besides the purity, which measures how mixed the noisy states become due to
losses, we shall quantify the disturbance of the pure cat states subjected to
the noise channel. For this, we use the fidelity between the pure and the
noisy states
$F(\outerproduct{\psi}{\psi},\hat{\rho}^{\eta})=\bra{\psi}\hat{\rho}^{\eta}\ket{\psi}$
to quantify how close the two states are jozsa1994fidelity . In Fig. 2 (b) we
show the fidelity between the noisy cat states with its pure counterpart for
the initial mean photon number of $\expectationvalue*{\hat{N}}=100$. Note that
for each $\eta$ we consider the pure state to be reduced in amplitude by a
factor of $\sqrt{\eta}$, in order to maximize the fidelity with the noisy
state which has a reduced amplitude due to the losses.

Figure 2: (a): Purity $\gamma^{\eta}$ of the noisy cat states
$\hat{\rho}^{\eta}$. (b): Fidelity between the noisy cat states and their pure
counterparts. The pure states have a reduced amplitude by a factor
$\sqrt{\eta}$ to maximize the fidelity. A vertical dashed line at $\eta=0.469$
(53.1% photon loss) indicates up to which $F_{\mathrm{HHG}}>F_{\pm}$. Note
that for $\eta\to 0$ all the noisy states become $\outerproduct{0}{0}$. In
both plots we consider an initial average photon number of
$\expectationvalue*{\hat{N}}=100$ before the noise channel, having $\alpha=10$
for the even/odd cat states, while $\alpha^{\prime}=10.5$ and
$\delta\alpha=-0.5$ for the HHG-cat.

The key results are that the HHG-cat state is significantly more robust than
the even/odd cat states for a wide range up to 53.1% photon loss. And we
emphasize, that the advantage of the HHG-cat state robustness is most
pronounced in the important regime of small photon losses, in which the
even/odd cat state counterparts exhibit significantly reduced purity and
fidelity and rapidly decay to the 50/50 mixture. In particular, for the
$\langle\hat{N}\rangle=100$ photon states considered in Fig. 2, the even/odd
cat state reduce to the maximally mixed state $\gamma_{\pm}^{\eta}=1/2$
already for $\eta=0.983$ while the HHG-cat nearly remains a pure state with
$\gamma^{\eta}_{\mathrm{HHG}}\approx 0.97$. In addition, the fidelity of the
even/odd cat state reduces to $F_{\pm}=1/2$ for $\eta=0.965$ while the noisy
HHG-cat is close to its pure counterpart with $F_{\mathrm{HHG}}\approx 0.97$.
In summary, the HHG-cat state is significantly more robust for high photon
numbers due to its unique structure. It is the relative displacement
$\delta\alpha$ in Eq. (1) which determines the distinguishability between the
two states in the HHG-cat, which remains constant even if the HHG-cat is
growing in size by means of the average photon number (due to increasing
$\alpha$). In contrast, for the even/odd cat states an increasing $\alpha$
would lead to an increased distinguishability, which cause a rapid decoherence
for high photon numbers already for very small losses of just a few photons.
This ultimately makes the even/odd cat states much more vulnerable compared to
the robust HHG-cat state. The HHG-cat state is protected against the losses
due to the indistinguishability of the two states in the superposition.

##  IV Metrological robustness

We shall now answer the question if the enhanced robustness of the HHG-cat
state against photon loss over the even/odd cat states carries over to the
problem of phase estimation leading to a metrological robustness against
decoherence. Like in the ideal case the QFI is used as the quantifier for the
metrological power of the cat states for the phase estimation problem. To
calculate the QFI for the noisy cat states, we use the general expression in
Eq. (2), and the eigenvalues from the spectral decomposition given in Eq. (6)
and Eq. (7). The QFI for the phase sensing using these noisy states is shown
in Fig. 3 for varying $\eta$, and normalized to the QFI for the pure odd cat
state for an average photon number of $\expectationvalue*{\hat{N}}=100$ before
losses. We can clearly see that the noisy HHG-cat state outperforms the noisy
odd cat state for a wide range of photon loss, and exhibit a bigger QFI up to
losses of 53.7% (indicated by a vertical dashed line). However, more
importantly, the main result of this work is that the noisy HHG-cat even
outperforms the noiseless pure odd cat state for up to 25.9% photon loss
(indicated by a vertical dashed line). These results imply that in a realistic
setting and for the specific phase detection setup we are considering, the
HHG-cat state provides an enhanced sensitivity compared to pure even/odd cat
states. Since experimentally relevant photon losses are usually much less than
25%, the HHG-cat state is significantly more useful than the noiseless
even/odd cat states, and most importantly, this even holds true when the HHG-
cat is subjected to losses. Finally, we observe that the noise range in which
the HHG-cat state is outperforming the even/odd cat states remains almost
constant over four orders of magnitude of the average photon number. This
shows that the metrological potential of the HHG-cat state is not restricted
to a particular number of photons in the field.

Figure 3: Quantum Fisher information for the noisy HHG-cat state, the pure and
noisy odd cat states, normalized to the QFI of the pure odd cat state
($\mathcal{F}/\mathcal{F}_{-}$) in the high photon number regime. Vertical
lines at $\eta=0.741$ (25.9% photon loss) and $\eta=0.463$ (53.7$\%$ photon
loss) indicate when $\mathcal{F}_{\mathrm{HHG}}^{\eta}\approx\mathcal{F}_{-}$
and $\mathcal{F}_{\mathrm{HHG}}^{\eta}\approx\mathcal{F}_{-}^{\eta}$,
respectively. The initial amplitude of the pure odd cat state is $\alpha=10$,
while for the HHG-cat we use $\alpha^{\prime}=10.5$ and $\delta\alpha=-0.5$
such that the average photon number is $\langle\hat{N}\rangle=100$. For the
noisy cases, the mean photon number for the LO is set to match the mean photon
number of the noisy cat states $|\beta|^{2}=100\eta$.

##  V Conclusions

We have demonstrated that the optical cat states obtained from quantum state
engineering approaches using the process of high harmonic generation
stammer2023quantum  show an enhanced robustness against photon losses compared
to the traditional even/odd cat states in the regime of high photon numbers.
Especially, the HHG-cat is significantly more robust in the regime of small
noise in which the loss of just a few photons from the initial 100 photons
brings the even/odd cat state into a maximally mixed state while the HHG-cat
almost remains a pure state. Together with the metrological advantage of the
HHG-cat in the ideal lossless scenario, this translates to a metrological
robustness for phase sensing measurements. In particular, we have shown that
the noisy HHG-cat state leads to a larger phase sensitivity than the noiseless
odd cat state even in the presence of more than $25\%$ photon loss, and
against the noisy odd cat state for up to more than $53\%$ losses. This is of
particular importance for protecting quantum technologies against noise
suter2016colloquium , and provides further insights into the applicability of
the connection between attosecond physics with quantum optics and quantum
information science bhattacharya2023strong ; lewenstein2022attosecond .

###### Acknowledgements.

P.S. acknowledges funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
847517. ICFO group acknowledges support from: ERC AdG NOQIA; MCIN/AEI
(PGC2018-0910.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan
National FIDEUA PID2019-106901GB-I00, Plan National STAMEENA
PID2022-139099NB-I00 project funded by MCIN/AEI/10.13039/501100011033 and by
the “European Union NextGenerationEU/PRTR” (PRTR-C17.I1), FPI); QUANTERA MAQS
PCI2019-111828-2); QUANTERA DYNAMITE PCI2022-132919 (QuantERA II Programme co-
funded by European Union’s Horizon 2020 program under Grant Agreement No
101017733), Ministry of Economic Affairs and Digital Transformation of the
Spanish Government through the QUANTUM ENIA project call – Quantum Spain
project, and by the European Union through the Recovery, Transformation, and
Resilience Plan – NextGenerationEU within the framework of the Digital Spain
2026 Agenda; Fundació Cellex; Fundació Mir-Puig; Generalitat de Catalunya
(European Social Fund FEDER and CERCA program, AGAUR Grant No. 2021 SGR 01452,
QuantumCAT U16-011424, co-funded by ERDF Operational Program of Catalonia
2014-2020); Barcelona Supercomputing Center MareNostrum (FI-2023-1-0013); EU
Quantum Flagship (PASQuanS2.1, 101113690); EU Horizon 2020 FET-OPEN OPTOlogic
(Grant No 899794); EU Horizon Europe Program (Grant Agreement 101080086 —
NeQST), ICFO Internal “QuantumGaudi” project; European Union’s Horizon 2020
program under the Marie Sklodowska-Curie grant agreement No 847648; “La Caixa”
Junior Leaders fellowships, La Caixa” Foundation (ID 100010434):
CF/BQ/PR23/11980043. Views and opinions expressed are, however, those of the
author(s) only and do not necessarily reflect those of the European Union,
European Commission, European Climate, Infrastructure and Environment
Executive Agency (CINEA), or any other granting authority. Neither the
European Union nor any granting authority can be held responsible for them.

## References

  * [1] V. Giovanneti, S. Lloyd, and L. Maccone.  Quantum-enhanced measurements: beating the standard quantum limit.  Science, 306:1330–1336, 2004.
  * [2] V. Giovannetti, S. Lloyd, and L. Maccone.  Advances in quantum metrology.  Nature Photonics, 5:222–229, 2011.
  * [3] J. Liu, H. Yuan, X.M. Lu, and X. Wang.  Quantum Fisher information matrix and multiparameter estimation.  Journal of Physics A: Mathematical and Theoretical, 53:023001, 2020\.
  * [4] C.W. Helstrom.  Quantum detection and estimation theory.  Journal of Statistical Physics, 1:231–252, 1969.
  * [5] G. Tóth and I. Apellaniz.  Quantum metrology from a quantum information science perspective.  Journal of Physics A: Mathematical and Theoretical, 47:424006, 2014\.
  * [6] J.J. Meyer.  Fisher information in noisy intermediate-scale quantum applications.  Quantum, 5:539, 2021.
  * [7] C.C. Gerry and P.L. Knight.  Quantum superpositions and Schrödinger cat states in quantum optics.  American Journal of Physics, 65:964–974, 1997.
  * [8] A. Gilchrist, K. Nemoto, W.J. Munro, T.C. Ralph, S. Glancy, S.L. Braunstein, and G.J. Milburn.  Schrödinger cats and their power for quantum information processing.  Journal of Optics B: Quantum and Semiclassical Optics, 6:S828, 2004\.
  * [9] T.C. Ralph, A. Gilchrist, G.J. Milburn, W.J. Munro, and S. Glancy.  Quantum computation with optical coherent states.  Physical Review A, 68:042319, 2003.
  * [10] J. Joo, W.J. Munro, and T.P. Spiller.  Quantum metrology with entangled coherent states.  Physical Review Letters, 107:083601, 2011.
  * [11] S.Y. Lee, C.W. Lee, H. Nha, and D. Kaszlikowski.  Quantum phase estimation using a multi-headed cat state.  Journal of the Optical Society of America B, 32:1186–1192, 2015\.
  * [12] M. Kitagawa and M. Ueda.  Squeezed spin states.  Physical Review A, 47:5138–5143, 1993.
  * [13] D.J. Wineland, J.J. Bollinger, W.M. Itano, and D.J. Heinzen.  Squeezed atomic states and projection noise in spectroscopy.  Physical Review A, 50:67–88, 1994.
  * [14] M. Płodzień, M. Kościelski, E. Witkowska, and A. Sinatra.  Producing and storing spin-squeezed states and Greenberger-Horne-Zeilinger states in a one-dimensional optical lattice.  Physical Review A, 102:013328, 2020.
  * [15] M. Płodzień, M. Lewenstein, E. Witkowska, and J. Chwedeńczuk.  One-axis twisting as a method of generating many-body Bell correlations.  Physical Review Letters, 129:250402, 2022.
  * [16] M. Płodzień, T. Wasak, E. Witkowska, M. Lewenstein, and J. Chwedeńczuk.  Generation of scalable many-body Bell correlations in spin chains with short-range two-body interactions.  arXiv:2306.06173, 2023.
  * [17] T. H. Yanes, M. Płodzień, M.M. Sinkevičienė, G. Žlabys, G. Juzeliūnas, and E. Witkowska.  One- and two-axis squeezing via laser coupling in an atomic Fermi-Hubbard model.  Physical Review Letters, 129:090403, 2022.
  * [18] M. Dziurawiec, T.H. Yanes, M. Płodzień, M. Gajda, M. Lewenstein, and E. Witkowska.  Accelerating many-body entanglement generation by dipolar interactions in the Bose-Hubbard model.  Physical Review A, 107, 2023.
  * [19] T.H. Yanes, G. Žlabys, M. Płodzień, D. Burba, M.M. Sinkevičienė, E. Witkowska, and G. Juzeliūnas.  Spin squeezing in open Heisenberg spin chains.  Physical Review B, 108:104301, 2023.
  * [20] B.M. Escher, R.L. de Matos Filho, and L. Davidovich.  General framework for estimating the ultimate precision limit in noisy quantum-enhanced metrology.  Nature Physics, 7:406–411, 2011.
  * [21] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier.  Generating optical Schrödinger kittens for quantum information processing.  Science, 312:83–86, 2006.
  * [22] A. Ourjoumtsev, H. Jeong, R. Tualle-Brouri, and P. Grangier.  Generation of optical ‘Schrödinger cats’ from photon number states.  Nature, 448:784–786, 2007.
  * [23] S. Glancy and H.M. de Vasconcelos.  Methods for producing optical coherent state superpositions.  Journal of the Optical Society of America B, 25:712–733, 2008.
  * [24] Y.M. Zhang, X.W. Li, W. Yang, and G.R. Jin.  Quantum fisher information of entangled coherent states in the presence of photon loss.  Physical Review A, 88(4):043832, 2013.
  * [25] X. Yu, X. Zhao, L. Shen, Y. Shao, J. Liu, and X. Wang.  Maximal quantum Fisher information for phase estimation without initial parity.  Optics Express, 26:16292–16302, 2018.
  * [26] M. Lewenstein, M.F. Ciappina, E. Pisanty, J. Rivera-Dean, P. Stammer, Th. Lamprou, and P. Tzallas.  Generation of optical Schrödinger cat states in intense laser–matter interactions.  Nature Physics, 17:1104–1108, 2021.
  * [27] P. Stammer, J. Rivera-Dean, A.S. Maxwell, Th. Lamprou, A.F. Ordóñez, M.F. Ciappina, P. Tzallas, and M. Lewenstein.  Quantum electrodynamics of intense laser-matter interactions: A tool for quantum state engineering.  Physical Review X Quantum, 4:010201, 2023.
  * [28] U. Bhattacharya, Th. Lamprou, A.S. Maxwell, A. Ordonez, E. Pisanty, J. Rivera-Dean, P. Stammer, M.F. Ciappina, M. Lewenstein, and P. Tzallas.  Strong–laser–field physics, non–classical light states and quantum information science.  Reports on Progress in Physics, 86:094401, 2023.
  * [29] J. Rivera-Dean, Th. Lamprou, E. Pisanty, P. Stammer, A.F. Ordóñez, A.S. Maxwell, M.F. Ciappina, M. Lewenstein, and P. Tzallas.  Strong laser fields and their power to generate controllable high-photon-number coherent-state superpositions.  Physical Review A, 105:033714, 2022.
  * [30] J. Rivera-Dean, P. Stammer, A.S. Maxwell, Th. Lamprou, P. Tzallas, M. Lewenstein, and M.F. Ciappina.  Light-matter entanglement after above-threshold ionization processes in atoms.  Physical Review A, 106:063705, 2022.
  * [31] P. Stammer, J. Rivera-Dean, Th. Lamprou, E. Pisanty, M.F. Ciappina, P. Tzallas, and M. Lewenstein.  High photon number entangled states and coherent state superposition from the extreme ultraviolet to the far infrared.  Physical Review Letters, 128:123603, 2022.
  * [32] P. Stammer.  Theory of entanglement and measurement in high-order harmonic generation.  Physical Review A, 106:L050402, 2022.
  * [33] P. Stammer, J. Rivera-Dean, A.S. Maxwell, Th. Lamprou, J. Argüello-Luengo, P. Tzallas, M.F. Ciappina, and M. Lewenstein.  Entanglement and squeezing of the optical field modes in high harmonic generation.  arXiv:2310.15030, 2023.
  * [34] J. Rivera-Dean, P. Stammer, E. Pisanty, Th. Lamprou, P. Tzallas, M. Lewenstein, and M.F. Ciappina.  New schemes for creating large optical Schrödinger cat states using strong laser fields.  Journal of Computational Electronics, 20:2111–2123, 2021.
  * [35] Th. Lamprou, J. Rivera-Dean, P. Stammer, M. Lewenstein, and P. Tzallas.  Nonlinear optics using intense optical Schrödinger “cat” states.  arXiv:2306.14480, 2023.
  * [36] P. Stammer.  On the role of the optical phase and coherence in high harmonic generation.  arXiv:2309.05010, 2023.
  * [37] J. Rivera-Dean, P. Stammer, A.S. Maxwell, Th. Lamprou, A.F. Ordóñez, E. Pisanty, P. Tzallas, M. Lewenstein, and M.F. Ciappina.  Quantum optical analysis of high-harmonic generation in solids within a Wannier-Bloch picture.  arXiv:2211.00033, 2022.
  * [38] J. Rivera-Dean, P. Stammer, A.S. Maxwell, Th. Lamprou, A.F. Ordóñez, E. Pisanty, P. Tzallas, M. Lewenstein, and M.F. Ciappina.  A Bloch-based quantum optical analysis of high-harmonic generation in semiconductors.  arXiv:2309.14435, 2023.
  * [39] A. Gorlach, M.E. Tzur, M. Birk, M. Krüger, N. Rivera, O. Cohen, and I. Kaminer.  High-harmonic generation driven by quantum light.  Nature Physics, pages 1–8, 2023.
  * [40] P. Stammer.  On the limitations of the semi-classical picture in high harmonic generation.  arXiv:2308.15087, 2023.
  * [41] A. Pizzi, A. Gorlach, N. Rivera, A. Nunnenkamp, and I. Kaminer.  Light emission from strongly driven many-body systems.  Nature Physics, 19(4):551–561, 2023.
  * [42] R.A. Campos, B.E.A. Saleh, and M.C. Teich.  Quantum-mechanical lossless beam splitter: SU(2) symmetry and photon statistics.  Physical Review A, 40:1371, 1989.
  * [43] M. Zhang, H. Kang, M. Wang, F. Xu, X. Su, and K. Peng.  Quantifying quantum coherence of optical cat states.  Photonics Research, 9:887–892, 2021.
  * [44] P.P. Rohde and A.P. Lund.  Practical effects in cat state breeding.  arXiv:quant-ph/0702064, 2007.
  * [45] A. Serafini, S. De Siena, F. Illuminati, and M.G.A. Paris.  Minimum decoherence cat-like states in Gaussian noisy channels.  Journal of Optics B: Quantum and Semiclassical Optics, 6:S591, 2004\.
  * [46] R. Jozsa.  Fidelity for mixed quantum states.  Journal of modern optics, 41:2315–2323, 1994.
  * [47] D. Suter and G.A. Álvarez.  Colloquium: Protecting quantum information against environmental noise.  Reviews of Modern Physics, 88:041001, 2016.
  * [48] M. Lewenstein, N. Baldelli, U. Bhattacharya, J. Biegert, M.F. Ciappina, U. Elu, T. Grass, P.T. Grochowski, A. Johnson, Th. Lamprou, et al.  Attosecond physics and quantum information science.  arXiv:2208.14769, 2022.