11institutetext: Universidade de São Paulo, IAG, Rua do Matão, 1226, Cidade
Universitária, 05508-090 São Paulo, Brazil
11email<EMAIL_ADDRESS>22institutetext: Rua Sessenta e Três,
125, Olinda, 53090-393 Pernambuco, Brazil
22email<EMAIL_ADDRESS>
# 3D stellar motion in the axisymmetric Galactic potential and the e–z
resonances

Tatiana A. Michtchenko 11 Douglas A. Barros 22

(Received —- –, —-; accepted —- –, —-)

###### Abstract

Context. The full phase-space information on the kinematics of a huge number
of stars provided by the Gaia Data Release 3 raises the demand for a better
understanding of the 3D stellar dynamics.

Aims. In this paper, we investigate the possible regimes of motion of stars in
the axisymmetric approximation of the Galactic potential, applying a 3D
observation-based model developed elsewhere. The model consists of three
components: the axisymmetric disk, the central spheroidal bulge and the
spherical halo of dark matter. The axisymmetric disk model is divided into
thin and thick stellar disks and H I and $\mathrm{H}_{2}$ gaseous disks
subcomponents, by combining three Miyamoto-Nagai disk proﬁles of any model
order (1, 2, or 3) for each disk subcomponent, to reproduce a radially
exponential mass distribution. The physical and structural parameters of the
Galaxy components are adjusted by observational kinematic constraints.

Methods. The phase space of the two-degrees-of-freedom model is studied by
means of the Poincaré and dynamical mapping, the dynamical spectrum method and
the direct numerical integrations of the Hamiltonian equations of motion.

Results. For the chosen physical parameters, the nearly-circular (close to the
rotation curve) and low-altitude stellar behaviour is composed of two weakly
coupled simple oscillations, radial and vertical motions. The amplitudes of
the vertical oscillations of these orbits are gradually increasing with the
growing Galactocentric distances, in concordance with the exponential mass
decay assumed. However, for increasing planar eccentricities $e$ and the
altitudes over the equatorial disk $z$, new regimes of stellar motion emerge
as a result of the beating between the radial and vertical oscillation
frequencies, which we refer to as $e$–$z$ resonances. The corresponding
resonant motion produces characteristic sudden increase/decrease of the
amplitude of the vertical oscillation, bifurcations in the dynamical spectra
and the chains of islands of stable motion in the phase-space.

Conclusions. The results obtained can be useful in the understanding and
interpretation of the features observed in the stellar 3D distribution around
the Sun.

###### Key Words.:

Galaxy: kinematics and dynamics – solar neighborhood – Galaxy: structure

##  1 Introduction

In the past decades, there has been a growing set of evidence that the disk
stars of the Milky Way galaxy exhibit density and velocity structures in
several phase-space planes of the Galactic disk, in both the radial and
vertical directions. The long-known stellar warp in the outer disk, a
vertically asymmetric distribution of stars close to the Galactic plane, is
seen as a bending of the plane upwards in the first and second Galactic
quadrants (longitudes $0^{\circ}\leq l\leq 180^{\circ}$) and downwards in the
third and fourth quadrants (longitudes $180^{\circ}\leq l\leq 360^{\circ}$)
(Momany et al., 2006). Recently, a number of stellar density substructures,
appearing as overdensities or dips, have been found in the LAMOST data, at
several radial and vertical positions in the Galactic disk (Wang et al.,
2018).

Large-scale vertical motions as a Galactic North–South asymmetry have been
revealed by spectroscopic and astrometric surveys such as SEGUE (Widrow et
al., 2012), RAVE (Williams et al., 2013), LAMOST (Carlin et al., 2013), as
well as the second Gaia data release (Gaia Collaboration et al., 2018). Such
bulk vertical motion of stars present a wave-like behaviour, which has been
referred to as bending and breathing motions, with stars coherently moving
towards or away from the Galactic mid-plane, on both sides of it (Kawata et
al., 2018; Ghosh et al., 2022; Khachaturyants et al., 2022). Several dynamical
processes have been evoked to explain these non-zero vertical motions: the
ones from external origins such as the passing of the Saggitarius dwarf galaxy
through the Milky Way disk (Gómez et al., 2013), or the excitation of the disk
due to interactions with dark matter subhaloes (Feldmann & Spolyar, 2015); and
the ones from non-axisymmetric internal perturbations such as the breathing
motion induced by the Galactic bar (Monari et al., 2015), or by the action of
the spiral density waves (Faure et al., 2014; Ghosh et al., 2022;
Khachaturyants et al., 2022, among others).

The density structures present in the $R$–$V_{\varphi}$ plane (Galactocentric
radius versus azimuthal velocity) of disk stars, known as diagonal ridges,
clearly visible in the distribution of the mean Galactic radial velocity
$\left<V_{R}\right>$, can also be seen in the vertical direction from maps of
the mean absolute distance from the disk mid-plane $\left<|z|\right>$ and the
mean vertical velocity $\left<V_{z}\right>$ of the stellar distribution
(Khanna et al., 2019; Wang et al., 2020). This attribute may be indicative of
a coupling between planar and vertical motions of the stars in the disk.
Several attempts to unravel the origin of these structures have been presented
in the literature, with some of them relying on simulations of external
mechanisms like the Sagittarius dwarf galaxy perturbation (Antoja et al.,
2018; Laporte et al., 2019; Khanna et al., 2019), and others from simulations
that take into account internal dynamics such as the bar or the spiral
resonances (Hunt et al., 2018; Michtchenko et al., 2018; Fragkoudi et al.,
2019; Barros et al., 2020).

Another striking feature associated with the vertical motion of the stars in
the Galactic disk is the phase-space spiral (or the so-called snail shell)
present in the $z$–$V_{z}$ plane. Many authors interpret it as evidence of an
ongoing phase mixing in the vertical direction of the disk due to an influence
of an outer perturbation (Antoja et al., 2018; Binney & Schönrich, 2018;
Bland-Hawthorn et al., 2019; Laporte et al., 2019), or caused by the
instability of a buckling bar (Khoperskov et al., 2019), or even due to the
earlier discussed vertical bending waves (Darling & Widrow, 2019) or vertical
breathing motion (Hunt et al., 2022, for two-armed phase spirals) in the disk.
Alternatively, Michtchenko et al. (2019) showed that the phase spiral can be
well-comprehended by the dynamical effects of the stellar moving groups in the
solar neighbourhood.

We suggest to analyse the above-cited Galactic structures and investigate
their causes by studying the stellar dynamics in the following steps:

  * •

Modeling a 3D Galactic gravitational potential. In this work, we use the
Barros et al. (2016) axisymmetric Galactic potential model, adjusted to recent
observational data. Firstly, the chosen potential is suitable to the study of
the vertical stellar motion and, secondly, it is simple in analytical forms
that gives a good representation of the radially exponential mass distribution
of the Galaxy. Moreover, any non-axisymmetric mass effects (e.g., due to the
central bar and/or spiral arms) can be easily added to the model posteriorly;

  * •

Studying the Galactic model in the whole phase space, in order to obtain all
the possible regimes of stellar motion that could provide reasonable
explanations to the observed Galactic phase-space structures. This is, in
part, the main objective of the present paper;

  * •

Performing numerical simulations through numerical integrations of stellar
orbits, to verify whether our achievements are satisfactory. This is a topic
for a future work.

In the present paper, we focus on the study of the stellar orbits in the
Galactic disk assuming a zeroth-order approximation of a time-independent and
axisymmetric Galactic gravitational potential. The whole phase space around
the Sun is investigated through techniques widely used in the Celestial
Mechanics. Resonances between the radial and vertical independent modes of the
stellar motion are characterized in terms of resonant zones, formed by islands
of stability that capture and trap stars inside them, enhancing the stellar
density, and separatrices that deplete the objects close to these regions. Our
goal, for a subsequent work, is to investigate possible associations between
the observed Galactic phase-space structures and the resonant zones originated
from the commensurabilities between the radial and vertical frequencies of the
stellar motion.

This paper is organized as follows. In Section 2, we describe the 3D
axisymmetric Galactic gravitational potential model used for the construction
of the Hamiltonian function and the calculus of the stellar orbits. In Section
3, we study the 3D stellar dynamics on the representative plane
$R$–$V_{\varphi}$ in the radial and vertical directions. In Section 4, regimes
of stellar motion are analysed in terms of the beating between the planar and
vertical frequencies of oscillation, and, in Section 5, we extend the analysis
for a wide range of initial vertical velocities $V_{z}$. Concluding remarks
are drawn in the closing Section 6.

##  2 3D Galactic model

Table 1: Physical and structural parameters of the disk components of the axisymmetric Galactic model Disk component | $M_{1}$ | $a_{1}$ | $M_{2}$ | $a_{2}$ | $M_{3}$ | $a_{3}$ | $b$
---|---|---|---|---|---|---|---
| ($10^{10}M_{\odot}$) | (kpc) | ($10^{10}M_{\odot}$) | (kpc) | ($10^{10}M_{\odot}$) | (kpc) | (kpc)
Thin disk | 2.282 | 3.859 | 2.342 | 9.052 | -1.846 | 3.107 | 0.243
Thick disk | 0.061 | 0.993 | 4.080 | 6.555 | -3.521 | 7.651 | 0.776
H I | 2.217 | 9.021 | 2.350 | 9.143 | -3.303 | 7.758 | 0.168
$\mathrm{H}_{2}$ | 1.005 | 6.062 | 0.177 | 3.141 | -0.907 | 4.485 | 0.128
Table 2: Physical and structural parameters of the spheroidal components of the Galactic model Bulge | $M_{b}$ | $a_{b}$
---|---|---
| ($10^{10}M_{\odot}$) | (kpc)
| 2.54 | 0.425
Dark halo | $r_{h}$ | $v_{h}$
| (kpc) | (km s-1)
| 5.56 | 169.77
Figure 1: Top: Rotation curve of the Galaxy. The observed rotation curve is
represented by the points with error bars, which indicate masers data from
high-mass star-forming regions (Reid et al., 2019; Rastorguev et al., 2017),
and H I and CO tangent-point data (Burton & Gordon, 1978; Clemens, 1985; Fich
et al., 1989). The red curve shows the analytical rotation curve expressed by
Eq. (10). The contributions of the modelled disk (continuous curve), bulge
(short-dashed line) and dark halo (long-dashed line) to the axisymmetric
Galactic potential are also shown. Bottom: Vertical force $K_{z}$ at $z=1.1$
kpc as a function of the Galactic radius. The black dots with error bars are
from observation measurements by Bovy & Rix (2013), while the grey solid curve
shows the $K_{z}$ force radial profile at $z=1.1$ kpc resulting from our
Galactic model.

Here, we briefly describe the 3D model for the axisymmetric Galactic potential
and refer the reader to Barros et al. (2016) for more details. The model
considers the contributions of three Galactic components into the potential;
they are the axisymmetric disk, the central spheroidal bulge and the spherical
halo of dark matter.

The axisymmetric disk is composed of the stellar (thin and thick disks) and
gaseous (H I and $\mathrm{H}_{2}$ disks) components. Each of these components
is a superposition of three Miyamoto–Nagai disks (MN, Miyamoto & Nagai, 1975),
which reproduces a radially exponential mass distribution (Smith et al.,
2015). We adopt the 1st–, 2nd– and 3rd–order expressions for the potential of
the MN–disks at the position $R$ and $z$ (cylindrical coordinates),
respectively:

$\Phi^{1}_{\mathrm{MN}}(R,z)=\frac{-GM}{\sqrt{R^{2}+(a+\zeta)^{2}}}\,,$ (1)
$\Phi^{2}_{\mathrm{MN}}(R,z)=\Phi^{1}_{\mathrm{MN}}(R,z)-\frac{GM\,a(a+\zeta)}{\left[R^{2}+\left(a+\zeta\right)^{2}\right]^{3/2}}\,,$
(2)
$\Phi^{3}_{\mathrm{MN}}(R,z)=\Phi^{2}_{\mathrm{MN}}(R,z)+\frac{GM}{3}\times\frac{a^{2}\left[R^{2}-2(a+\zeta)^{2}\right]}{\left[R^{2}+\left(a+\zeta\right)^{2}\right]^{5/2}}\,,$
(3)

where $\zeta=\sqrt{z^{2}+b^{2}}$, with $b$ being the vertical scale length,
while $M$ and $a$ are the mass and the radial scale length of each disk
component, respectively.

The potential of the thin disk is then written as a combination of three
third-order MN–disks

$\Phi_{\mathrm{thin}}(R,z)=\sum_{i=1}^{3}{\Phi^{3}_{\mathrm{MN}}}_{i}(R,z)\,,$
(4)

while the potential of the thick disk is a combination of three first-order
MN–disks

$\Phi_{\mathrm{thick}}(R,z)=\sum_{i=1}^{3}{\Phi^{1}_{\mathrm{MN}}}_{i}(R,z)\,.$
(5)

The contributions of the gaseous disks components H I and $\mathrm{H}_{2}$
into the total axisymmetric potential are, respectively,

$\displaystyle\Phi_{\mathrm{HI}}(R,z)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{3}{\Phi^{2}_{\mathrm{MN}}}_{i}(R,z)\,,$ (6)
$\displaystyle\Phi_{\mathrm{H_{2}}}(R,z)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{3}{\Phi^{3}_{\mathrm{MN}}}_{i}(R,z)\,.$ (7)

The potential of the Galactic bulge is derived by a Hernquist density
distribution proﬁle, in the form (Hernquist, 1990):

$\Phi_{\mathrm{b}}(R,z)=\frac{-G\,M_{\mathrm{b}}}{\sqrt{R^{2}+z^{2}}+a_{\mathrm{b}}}\,,$
(8)

where two free parameters, $M_{\mathrm{b}}$ and $a_{\mathrm{b}}$, are the
total mass and the scale radius of the bulge, respectively.

Finally, we consider a spherical dark halo, whose potential is modelled with a
logarithmic potential in the form (e.g. Binney & Tremaine, 2008):

$\Phi_{\mathrm{h}}(R,z)=\frac{v_{\mathrm{h}}^{2}}{2}\,\ln\left(R^{2}+z^{2}+r_{\mathrm{h}}^{2}\right)\,,$
(9)

where $r_{\mathrm{h}}$ is the core radius and $v_{\mathrm{h}}$ is the circular
velocity at large $R$ (i.e., relative to the core radius), respectively.

The analytical rotation curve is then calculated using the expression.

$V_{rot}(R)=\sqrt{R\times\frac{\partial\Phi_{0}}{\partial R}|_{z=0}}\,,$ (10)

where the axisymmetric potential $\Phi_{0}$ is just the sum of the potentials
of the Galactic components described above, that is:

$\Phi_{0}(R,z)=\Phi_{\mathrm{thin}}+\Phi_{\mathrm{thick}}+\Phi_{\mathrm{HI}}+\Phi_{\mathrm{H_{2}}}+\Phi_{\mathrm{b}}+\Phi_{\mathrm{h}}\,.$
(11)

The physical and structural parameters of the components of the Galactic
axisymmetric potential are obtained by applying a fitting procedure, which
adjusts the analytical rotation curve Eq. (10) to the observed one. For the
observation-based rotation curve, we use data of H I-line tangential
directions from Burton & Gordon (1978) and Fich et al. (1989), CO-line
tangential directions from Clemens (1985), and maser sources data associated
with high-mass star-forming regions from Reid et al. (2019) and Rastorguev et
al. (2017). From these data, the Galactic radii and rotation velocities were
calculated taking the distance of the Sun from the Galactic center as
$R_{0}=8.122$ kpc (GRAVITY Collaboration et al., 2018) and the Local Standard
of Rest (LSR) velocity at the Sun $V_{0}=233.4$ km s-1 (Drimmel & Poggio,
2018). As additional constraints to the fitting procedure, we used the value
of the local angular velocity $\Omega_{0}$, the local disk surface density
$\Sigma_{0}$, and the local surface density within $|z|\leq 1.1$ kpc. For
details of the data used and the fitting procedure, see Barros et al. (2016,
2021).

The parameters obtained are shown in Table 1 and Table 2. The contributions of
each component of the axisymmetric potential to the Galactic rotation curve
are plotted in Fig. 1 (top panel), together with the calculated rotation curve
(red line) and the data points. Figure 1 (bottom panel) also shows the radial
profile of the $K_{z}$ vertical force resulting from the Galactic model
calculated at $z=1.1$ kpc. We note a good agreement between the modeled force
and the measured $K_{z}$–force, based on measurements by Bovy & Rix (2013)
(black dots with error bars).

The stellar orbits are calculated through the numerical integrations of the
equations of motions defined by the Hamiltonian function given as

$H(R,p_{r},z,V_{z})=\frac{1}{2}\left[p_{r}^{2}+\frac{L_{z}^{2}}{R^{2}}+V_{z}^{2}\right]+\Phi_{0}(R,z)\,,$
(12)

where $p_{r}$ and $V_{z}$ are the linear radial momentum and vertical
velocity, respectively, while $L_{z}$ is the angular momentum, which is a
constant in the axisymmetric approximation (all momenta are given per mass
unit)

##  3  3D stellar dynamics on the representative plane

Figure 2: Topology of the Hamiltonian (Eq. 12) on the representative
$R$–$V_{\varphi}$ plane. Continuous curves are the levels of the angular
momentum $L_{z}=R\times V_{\varphi}$, while dashed curves are the energy
levels, calculated with the fixed initial values $p_{R}=0$, $z=0$ and
$V_{z}=20$ km s-1. Rotation curve is shown by red dots. The projections of the
three 3D orbits on the plane are: the Sun’s orbit (cyan dots), one orbit
starting at the configuration 1 and other at the configuration 2 (black dots).
Note that the motion occurs always along a fixed $L_{z}$-value, due to
conservation of the angular momentum $L_{z}$. The conservation of the energy
during the motion defines two turning points of each orbit, both lying at the
intersections of the corresponding $L_{z}$– and energy levels. Top panel: the
total energy (see Eq. 12), in arbitrary units, calculated along the rotation
curve, as a function of $R$.

The dynamical model defined by the Hamiltonian in Eq. (12) is of two-degree-
of-freedom. The resulting stellar motion can be formally represented by two
coupled oscillations: one, in the radial (equatorial) direction, is described
by the pair of the variables $R$–$p_{r}$, and other, in the vertical
direction, described by the pair $z$–$V_{z}$. The phase space of the system is
four-dimensional, that makes it difficult to visualize the stellar orbits.
Yet, in this paper, we introduce a representative plane, which allows us to
present the main features of the 3D stellar dynamics. This plane is the
$R$–$V_{\varphi}$ plane, where $V_{\varphi}$ is the tangential velocity
defined as $V_{\varphi}=L_{z}/R$; we show this plane in Fig. 2.

To study the stellar motions on the $R$–$V_{\varphi}$ plane, we fix, in this
paper, the initial values of the planar linear momentum at $p_{r}=0$ and the
vertical height at $z=0$. It is worth noting that this choice preserves the
generality of the presentation, since all closed stellar motions under the
potential given by Eq. (11) pass through these conditions. The initial value
of the vertical velocity is chosen as $V_{z}=20$ km s-1, except when the
dependence of the motion on the initial vertical configuration is analysed
(Sect. 5). The chosen velocity value is close to the value of the standard
deviation of the $V_{z}$ distribution of the stars from the Gaia DR3 (Gaia
Collaboration et al., 2016, 2023), with $|z|<1$ pc (we obtained
$\sigma_{V_{z}}=18$ km s-1).

First, we plot on the $R$–$V_{\varphi}$ plane in Fig. 2 the main
characteristics of the Hamiltonian dynamics: the orbital energy and the
angular momentum, which are both constants of motion in our model. The energy
levels and $L_{z}$–levels are shown by dashed and continuous lines,
respectively. The rotation curve obtained using Eq. (10), is shown by the red
curve.

By definition, the rotation curve is a place of the circular orbits, which are
orbits of the minimal energy, for a given $L_{z}$-value. The geometrical
interpretation of this condition is that, given the location of a circular
orbit on the rotation curve, the corresponding $L_{z}$-level is tangent to the
minimal energy level on the $R$–$V_{\varphi}$ plane. This fact can be observed
in Fig. 2, considering that the energy is continuously increasing with the
increasing radial distances $R$, as shown on the top panel in Fig. 2. The same
$L_{z}$-level intersects the levels of higher energy always at two points,
which are turning points of the radial oscillations at the condition
$p_{r}=0$.

We plot in Fig. 2 the direct projections of three 3D stellar orbits: the one
is the Sun’s orbits (cyan dots) calculated with initial values $R=8.122$ kpc,
$p_{r}=-12.9$ km s-1, $V_{\varphi}=245.6$ km s-1, $z=0.015$ kpc and
$V_{z}=7.78$ km s-1 (Drimmel & Poggio, 2018). The orbits of other two
fictitious stars (black dots) were starting at points 1 and 2, respectively,
with $p_{r}=0$, $z=0$ and $V_{z}=20$ km s-1. We can verify that all orbits are
closed, that is, each one evolves between two turning points, which belong to
the same energy levels, along the corresponding $L_{z}$–levels (constant of
motion), around the circular orbits, which lie at the intersection of the
corresponding energy and $L_{z}$–levels.

The detailed analysis of the main features of the stellar motions on the
representative $R$–$V_{\varphi}$ plane can be done in two steps: first,
focusing on the radial oscillations on the Galactic equatorial plane, and,
then, on the vertical oscillations around the equatorial plane. By
understanding the behaviour of each component of the stellar motion in the
axisymmetric potential, we can consider new effects produced by their
interaction. We will do it in the next sections.

Figure 3: Poincaré map of the stellar 3D orbits calculated with the initial
conditions chosen along the Sun’s $L_{z}$–level. The initial vertical
conditions were chosen at $z=0$ and $V_{z}=20$ km s-1. The projections are
calculated at instants when the star crosses the equatorial plane with
$V_{z}>0$. The orbits are concentric, with the circular orbit at the center
(black dot). This orbit lies at the intersection of the rotation curve and the
Sun’s $L_{z}$–level on the $R$–$V_{\varphi}$ plane, at $R_{c}=8.522$ kpc and
$V_{\varphi}=234$ km s-1. Red and cyan dots show the projections of two
periodic orbits on the Poincaré map corresponding to the 2/1 and 4/1
resonances, respectively (see later Fig. 7).

###  3.1 Radial oscillations

The $R$–$V_{\varphi}$ plane shown in Fig. 2 is commonly used for understanding
the radial (or planar) motion of stars. Indeed, due to the conservation of the
angular momentum during the motion, the projections of the 3D orbits are
aligned along the corresponding $L_{z}$–levels (continuous lines), as observed
in the case of our three examples. Since the orbital energy is also conserved
in our problem, the orbit’s projections must match the corresponding energy
level at conditions $p_{r}=0$, which was chosen in the construction of the
plane. By definition, this condition defines two turning points of a closed
orbit, which delimit the path of the radial oscillation on the
$R$–$V_{\varphi}$ plane and, consequently, the maximal ($R_{\mathrm{max}}$)
and minimal ($R_{\mathrm{min}}$) values of the radial orbital variation (and
the tangential velocity $V_{\varphi}$, for a given $L_{z}$-value).

Any $L_{z}$–level crosses the rotation curve (red curve) at one point at
$R_{c}$ on the $R$–$V_{\varphi}$ plane in Fig. 2, which corresponds to the
circular orbit of the minimal energy for the corresponding $L_{z}$. At this
condition, two turning points merge at one point, characterizing the circular
orbit with the zero-amplitude $R$–oscillation, that is, the periodic orbit of
the two-degree-of-freedom problem. All other orbits calculated along the same
$L_{z}$–level, are quasi-periodic orbits oscillating with the non-zero
$R$–amplitudes, in such a way that larger the deviations from the rotation
curve, larger are the amplitudes of the radial variations.

The amplitude of the radial oscillation is related to the planar eccentricity
of the orbit as

$e=\frac{R_{\mathrm{max}}-R_{\mathrm{min}}}{2R_{c}}\,,$

which is uniquely defined by the values of the orbital energy and angular
momentum of the stellar orbit. In the conservative problem considered in this
paper, it is a constant of motion. The gain/loss of the orbital energy due to
some external physical processes, without changing the angular momentum of the
system, will increase/decrease the amplitude of $R$-oscillation and,
consequently, the planar eccentricity. This effect, known as ‘heating’ in the
stellar dynamics, is analogous to the tidal gravitational interactions in the
planet-satellite systems.

The equatorial motions of the stars are shown in Fig.3, where we plot the
orbits in the subspace $R$–$p_{R}$. The orbits were calculated through
numerical integrations of the equations of motion defined by the Hamiltonian
(Eq. 12), with initial conditions chosen along the level $L_{z}=1995$ kpc km
s-1, which corresponds to the Sun’s orbit (see Fig. 2). To obtain the
projection of the 3D motion on the Galactic plane, we gather the orbital
coordinates at the instants when the orbits cross the equatorial plane (the
condition $z=0$), in the direction of the positive $z$-values (the condition
$V_{z}>0$). This approach is known as the Poincaré surfaces of section method,
which allows us to plot the radial motion separately from the vertical one.

Figure 4: Representative $R$–$V_{\varphi}$ plane shown in Fig. 2, with the
same rotation curve (red dots) and the Sun’s orbit (cyan dots) along the
$L_{z}$–level (dashed line). The equidistant levels are used to represent the
amplitude of the vertical oscillation of the stellar orbits (see top panel to
associate the absolute $z_{\mathrm{max}}$-values), with the initial conditions
$z=0$ and $V_{z}=20$ km s-1. Top panel: the amplitude of the vertical
oscillation, $z_{\mathrm{max}}$, calculated along the rotation curve, as a
function of $R$.

Figure 3 shows that all orbits in the $R$–$p_{R}$ subspace oscillate around a
circular periodic orbit (black dot), which lies at the intersection of the
corresponding $L_{z}$–level and the rotation curve, at $R_{c}=8.522$ kpc and
$V_{\varphi}=234$ km s-1 on the $R$–$V_{\varphi}$ plane in Fig. 2. This
periodic orbit appears as a fixed point on the $R$–$p_{R}$ plane, indicating
that the amplitude of the radial oscillation is zero (or $e=0$). At the same
time, the vertical motion started with the initial $V_{z}=20$ km s-1 has non-
zero amplitude; it is a simple oscillation around the equatorial plane
($z=0$).

In the vicinity of the fixed point (black dot), between 7.0 kpc and 10.5 kpc,
the concentric low-eccentricity orbits, with $e<0.2$, evolve in good agreement
with the epicyclic approximation (e.g. Binney & Tremaine, 2008). The orbit of
the Sun, with eccentricity of 0.061, is located in this region of the phase
space. The values of the radial and vertical frequencies for the Sun’s $L_{z}$
level, derived under the epicyclic approximation, are equal to 0.00542 and
0.01387 (in units of 1/Myr), respectively, that corresponds to the radial and
vertical periods of 184.5 Myr and 72.07 Myr. These periods derived from the
power spectrum of the Sun’s orbit, calculated numerically through our model,
are of 161.3 Myr and 75.2 Myr, respectively.

The orbits with the increasing eccentricity, however, deviate from the
epicyclic approximation and, when the eccentricity approaches 0.5, a notable
feature appears in Fig.3. Indeed, at the radial distances beyond 12 kpc, we
can clearly observe a lack of circulating orbits, but there are some islands
on the Poincaré map. As shown below, the orbits in this region are different
structurally from the low-eccentricity orbits and are related to the new kind
of motion mode. The planar axisymmetric model is not able to explain this
feature, thus, we proceed investigating the vertical motion and its
interaction with the planar component.

Figure 5: Same as in Fig. 3, except the projections were calculated at
conditions $p_{R}=0$ and $\dot{p}_{R}>0$ on the $z$–$V_{z}$ plane. Red, cyan
and green dots show the projections of three periodic orbits on the Poincaré
map corresponding to the 2/1, 4/1 and 6/1 resonances, respectively. The large
island of the orbits oscillating around the center at $z=0$ kpc and
$V_{z}=103.0$ km s-1 (black dot) is related to the strong 1/1 resonance (see
later Fig. 7).

###  3.2 The vertical oscillations

The same representative $R$–$V_{\varphi}$ plane shown in Fig. 2 can be used to
study the stellar motion in the vertical direction. For this, we plot in Fig.
4, by continuous lines, the levels of the maximal altitude over the equatorial
plane that a star reaches during its vertical oscillation
($z_{\mathrm{max}}$). To obtain $z_{\mathrm{max}}$, we integrate numerically,
over several billions years, the equations of motion defined by the
Hamiltonian (Eq. 12), over the 200$\times$200 grid of the initial conditions
covering the $R$–$V_{\varphi}$ plane. All orbits started with $p_{R}=0$ on the
galactic equatorial plane (at $z=0$) and with the same initial vertical
velocity $V_{z}=20$ km s-1. The panel on the top of Fig. 4 shows
$z_{\mathrm{max}}$ obtained in this way for the circular orbits located along
the rotation curve (red curve). We plot also the projection of the Sun’s orbit
(cyan dots) along the $L_{z}$–level (dashed line) on the $R$–$V_{\varphi}$
plane.

To understand the $z_{\mathrm{max}}$–evolution on the Galactic
$R$–$V_{\varphi}$–plane, we consider first the nearly circular orbits,
distributed along the rotation curve. For these low-eccentricity orbits, the
amplitude of the vertical oscillation is smoothly increasing with the growing
Galactocentric distance, as shown on the top panel in Fig. 4. This result is
consistent with the radially decaying exponential mass distribution model
applied and is in agreements with the observed increase of the scale height
with the increase of the Galactic radius of the disk stars, which is the well-
known flaring of the Galactic disk (López-Corredoira et al., 2002; Amôres et
al., 2017; Robin et al., 2022).

However, the continuous evolution of the $z_{\mathrm{max}}$–levels is
interrupted in the regions of the increasing planar eccentricities when we
move away from the rotation curve on the $R$–$V_{\varphi}$ plane in Fig. 4, in
the direction of both the higher and lower tangential velocities. To
understand this behaviour, it is worth noting that, in domains of regular
oscillations, small changes in the initial conditions lead to small changes of
the elements of the orbits, in particular, the amplitude $z_{\mathrm{max}}$.
Consequently, the levels suffer slight displacements on the map when the
initial conditions are gradually changed. However, in the vicinity of the
resonances, small changes in the initial configurations produce large changes
of the orbital elements, forming singular structures, such as resonant islands
and stochastic layers. On the dynamical map in Fig. 4, these structures appear
in the form of stalactites of different widths.

Figure 6: Top: Amplitude of the vertical oscillation of stars as a function of
the initial Galactocentric distance. The initial values of $V_{\varphi}$ are
chosen along the Sun’s $L_{z}$–level, while $z=0$ and $V_{Z}=20$ km s-1.
Bottom: Dynamical spectrum of the stellar orbits calculated along the Sun’s
$L_{z}$–level. The nominal positions of some resonances are shown by vertical
dashed lines.

On the Poincaré map in Fig. 5, the same structures appear as chains of islands
contrasting with the orbits, which regularly circulate around the origin. To
construct this map, we pick up the orbital coordinates $z$ and $V_{z}$ at the
instants when the star is in the turning point of its orbit ($p_{r}=0$) of the
minimal radial distance ($\dot{p}_{r}>0$). The initial configurations of the
orbits were chosen along the Sun’s $L_{z}$–level and $z$ and $V_{z}$ fixed at
0 and 20 km s-1, respectively, as in Fig. 3. However, in this case, the
initial $R$–values were extended to the very large radial distances, up to 27
kpc.

Comparing this map to the map of the equatorial motion shown in Fig. 3, we
verify that the behaviour in the vertical direction is more complicated. We
can observe the projections of the qualitatively different types: there are
circulating orbits of regular motion at smaller stellar altitudes. With the
growing height, the chains of islands of different width appear; they are
separated by the stochastic layers. A large island dominates the region of
high vertical velocities in the positive half-plane of the map and is
surrounded by the sea of chaotic motion. To identify the cause of this
behaviour, we apply a special spectral analysis method and describe the
results obtained in the next section.

Figure 7: Resonant periodic orbits on the $R$–$p_{R}$ plane (top row) and
$z$–$V_{z}$ plane (bottom row). The initial conditions were chosen along the
Sun’s $L_{z}$–level, while $z=0$ and $V_{z}=20$ km s-1. Column (a): the 2/1
resonant orbit with $e=0.45$ (red dots on the Poincaré maps in Figs. 3 and 5).
Column (b): the 4/1 resonant orbit with $e=0.72$ (cyan dots in Figs. 3 and 5).
Column (c): the 6/1 resonant orbit with $e=1.01$ (green dots in Fig. 5).
Column (d): the 1/1 resonant orbit with $e=0.802$ (black dot in Fig. 5).

##  4 The e–z resonances

To understand the features of the resonant motion, we start plotting the
amplitude of the vertical oscillation ($z_{\mathrm{max}}$) as a function of
the initial radial distance $R$ on the top panel in Fig. 6. In contrast to one
another shown on the top panel in Fig. 4, the amplitude $z_{\mathrm{max}}$ was
now calculated with initial conditions chosen along the $L_{z}$–level
corresponding to the Sun’s orbit. In this case, all orbits are eccentric, with
exception of one located at $R_{c}=8.522$ kpc, at the intersection of the
$L_{z}$–level with the rotation curve in Fig. 4.

For the orbits starting along the same $L_{z}$–level at the same initial
configurations in the vertical direction ($z=0$ and $V_{z}=20$ km s-1), the
minimal value of $z_{\mathrm{max}}$ is associated with the circular orbit, at
$R_{c}=8.522$ kpc. The vertical amplitude is generally increasing when the
orbital eccentricity grows with both increasing and decreasing distances from
the circular orbit. However, this evolution of $z_{\mathrm{max}}$ is not
monotonous, but shows sudden increase/decrease at some radial distances. The
nature of this behaviour can be analyzed using the dynamical power spectrum
method, which allows us to detect the change of a regime of motion through the
analysis of the evolution of the main frequencies of the dynamical system
(detailed description of the method and its applications to different systems
are found in e.g. Michtchenko et al., 2002, 2017).

Figure 6 bottom shows two proper frequencies, $f_{r}$ and $f_{z}$, which are
the frequencies of the radial and vertical oscillations, respectively, as
functions of the radial distances $R$. The evolution of the $f_{r}$–frequency
(red dots) and its harmonics is monotonous over the whole $R$–range; it
reaches the maximal value at $R_{c}=8.522$ kpc, that is expected since the
circular orbit is an equilibrium solution of the effective potential
$\Phi_{\mathrm{eff}}=\frac{L_{z}^{2}}{2R^{2}}+\Phi_{0}(R,z)$, with the
constant $L_{z}$. On the other hand, the behaviour of the $f_{z}$–frequency
(black dots) on the dynamical spectrum is highly non-harmonic and we can
observe the appearance of bifurcations, stable islands and an erratic
scattering of the dots, which is characteristic of the chaotic motion.

Analyzing the evolution of the proper frequencies in Fig. 6 bottom, we verify
that bifurcations occur when the beats between the two independent frequencies
occur. This condition is known as resonances between distinct modes of motion,
which lead to changes in the regime of motion. One beat occurs around $R=14$
kpc, where $f_{R}\cong 2\,f_{z}$, giving origin to the island of the 2/1
resonant motion. The periodic orbit of the 2/1 resonance is shown in the
subspaces $R$–$p_{R}$ (top) and $z$–$V_{z}$ (bottom) in Fig. 7 (column a); the
projections of this orbit in the Poincaré maps are shown by red dots in Figs.
3 and 5.

The beating frequencies are also observed at $R=17.8$ kpc in the dynamical
spectrum in Fig. 6 bottom; this event is associated to the 4/1 resonance. The
radial and vertical oscillation modes of the corresponding periodic orbit are
shown in Fig. 7 (column b); the projection of this orbit in the Poincaré maps
is shown by cyan dots in Figs. 3 and 5. In Fig. 5, we also show, by green
dots, the projection of the 6/1 resonant orbit (see Fig. 7 c); the location of
this resonance in the dynamical spectrum is shown at $R=19.2$ kpc.

We denote these resonances as the $e$–$z$ resonances, where $e$ and $z$ denote
planar eccentricities and vertical heights of stellar orbits, respectively.
The most important from the $e$–$z$ resonances is the 1/1 resonance, whose
domain is large and separated by the layers of chaotic motion (separatrices).
For a given value of $L_{z}$, it appears at very large Galactic distances; in
the dynamical map in Fig. 6, its domain is extended from 20 kpc to 35 kpc. The
motion inside the 1/1 resonance is stable, with the amplitude of the vertical
oscillation reaching its minimal value at the center (locus) of the resonance,
at $\sim$27 kpc; this locus appears as a fixed point (black dot) in the
Poincaré map in Fig. 5. The locus is a periodic orbit of the 1/1 resonance
shown in Fig. 7 (column d). The neighborhood of the 1/1 resonance is generally
a domain of highly unstable motion. This is due to the overlap of the high-
order resonances of the type $f_{r}/f_{z}\cong m/n$ ($m$ and $n$ are
integers).

Finally, it is worth emphasizing that, comparing the evolution of the two
frequencies in Fig. 6 bottom, we note the qualitative difference in the
behaviour of two independent modes: while the radial motion seems to be
unaffected by the passages through the resonances, their impact on the
vertical motion is significant. This could be explained by the peculiar
characteristics of the massive disk potential.

##  5 The dependence on the initial vertical velocity $V_{z}$

Figure 8: The $R$–$V_{z}$ plane of initial conditions chosen along the Sun’s
$L_{z}$-level, with $p_{R}=0$ and $z=0$. The equidistant levels of the maximal
vertical deviation from the equatorial plane, $z_{\textrm{max}}$, are shown by
continuous curves; its value is mainly increasing with the increasing $V_{z}$
(see beside panel to associate the absolute $z_{\mathrm{max}}$-values). The
red curve shows the location of the circular orbits with $L_{z}=1995$ kpc km
s-1 as a function of the initial $V_{z}$. The structures, which are
characteristic of the $e$–$z$ resonances, appears on both sides of the red
curve, expanding their domains with increasing eccentricities. The blue curve
shows loci of the very strong 1/1 resonance. The projection of the Sun’s orbit
is shown by cyan dots; its maximal $z$–amplitude is around 0.85 kpc. Beside
panel: $z_{\mathrm{max}}$-values as a function of the initial $V_{z}$, for
circular orbits (red curve) and orbits with the initial $R=13$ kpc (black
curve).

In this section, we investigate the resonant behaviour as a function of the
initial vertical velocity, $V_{z}$. For this, we introduce the representative
plane $R$–$V_{z}$ of the initial conditions and fix the rest of the variables
at $p_{r}=0$, $z=0$ and $L_{z}=1995$ kpc km s-1, this last corresponding to
the Sun’s angular momentum.

Figure 8 presents the amplitude of the vertical oscillation in the form of the
$z_{\textrm{max}}$–levels on the representative $R$–$V_{z}$ plane. The red
curve shows the positions of the circular orbits of the constant $L_{z}$,
which are slightly dislocated in the direction of the larger Galactocentric
distances, when the initial vertical velocities increase. The panel beside the
$R$–$V_{z}$ plane allows us to quantify the $z_{\textrm{max}}$–amplitude,
showing its values calculated along the loci of the circular orbits (red
curve). The smoothed evolution of this quantity with the increasing initial
$V_{z}$ can be observed.

The bifurcation structures, which are characteristic of the resonances, appear
outside the loci of the circular orbits in Fig. 8 and are strengthened with
the increasing planar eccentricities of the orbits. The black curve on the
beside panel shows the evolution of the $z_{\textrm{max}}$–amplitude
calculated along the constant initial $R=13$ kpc. We can observe the behaviour
which is similar to that shown on the top panel in Fig. 6, indicating the
passages through some $e$–$z$ resonances. The most prominent passage is
through the 1/1 resonance, whose loci are shown by the blue dots in Fig. 8.
The projection of the Sun’s trajectory on the $R$–$V_{z}$ plane is shown by
the cyan dots. Its proximity to the rotation curve avoids the capture in one
of the resonances.

##  6 Conclusions

In this paper, we report the existence of the resonant motion of the kind
$e$–$z$ ($e$ and $z$ being the planar eccentricity and the vertical altitude
of stellar orbits, respectively), produced by the axisymmetric Galactic 3D
potential.

We investigate the spacial motion of the stars applying the elaborated model
for the axisymmetric potential of the Galactic disk (Barros et al., 2016). The
model accounts for the combined gravitational effects due to the two stellar
disks (thin and thick disks) and two gaseous disks (H I and $\mathrm{H}_{2}$
disks). Moreover, to account for a radially exponential Galactic mass
distribution, each of the disks is approximated by the composition of the
three Miyamoto-Nagai disk proﬁles (Smith et al., 2015). The physical and
structural parameters of the modeled components were chosen to correspond to
the observable rotation curve of the Galaxy.

The motion of stars in the axisymmetric Galactic potential is investigated in
the whole phase space applying several techniques from the Celestial
Mechanics. The one of those is the dynamical mapping of the representative
plane chosen here as the $R$–$V_{\varphi}$ plane. Based on the conservation
laws, we identify two independent modes of the stellar motion, radial and
vertical ones. The radial motion is an oscillation around the circular
solution, which is a circular orbit belonging to the rotation curve and
characterized by the same angular momentum of the radial mode. The vertical
motion is an oscillation around the equatorial Galactic plane.

For nearly circular orbits, with small planar eccentricities, two oscillations
are weakly coupled. It is worth noting that the circular orbits survive even
in the case when vertical velocities are very high. However, when the planar
eccentricities increase, the coupling between two modes of motion also
increase. Their interaction is better visualised in the dynamical power
spectra, which clearly show the regions in the phase space where two proper
frequencies are commensurable. This beating frequencies indicate the resonant
zones, which we denote as $e$–$z$ resonances.

In the Hamiltonian theories, resonances are fundamental properties of
dynamical systems with two or more degrees of freedom. They occur in domains
of the phase space where the frequencies of the independent modes of motion
are commensurable. The locations and the sizes of the resonant domains are
strongly dependent on the physical parameters of the model adopted to describe
the system under study. Thus, associating some observable structures of the
objects to the resonances, we can assess the reliable ranges of the parameters
of the applied model. One example of this approach is shown in Michtchenko et
al. (2018), where we relate the moving groups in the solar neighbourhood to
the Lindblad resonances produced by the spirals perturbations, and estimate
the spiral strength and pattern speed.

The main features of the $e$–$z$ resonant motion are sudden increase/decrease
in the amplitude of the vertical oscillation, capture and trap of the stars
inside the stable resonant zones, enhancing the density, and deplete of the
regions close to the separatrices of the resonances. This behaviour is
characteristic of the vertical oscillations, while the planar radial motion is
only slightly affected by the $e$–$z$ resonances; this is probably due to
peculiar properties of a disk potential at all.

It is worth emphasizing that, despite the fact that the resonances occur at
very large Galactic distances, the high eccentricities of the resonant orbits
allow us their detection in the solar neighbourhood, according to Fig. 7.
Indeed, analyzing the behaviour of the Gaia-eDR3 sample from the solar
vicinity ($R_{\odot}\pm 0.2$ kpc and $|p_{R}|\geq 100$ km s-1), we have
detected 10% (from one thousand randomly chosen stars) of the objects showing
the resonant oscillations, mainly, inside the 2/1 resonance. Of course,
whether these stars are really evolving inside the $e$–$z$ resonances, it
depends strongly on the parameters adopted in our model. Thus, we should be
able to observe the described manifestations of the $e$–$z$ resonances
analysing the distributions of the observable proper elements of the stars.
The comparison to theoretical predictions could allow us to improve unknown
values of the parameters, which describe the observable Galactic mass
distribution.

###### Acknowledgements.

We acknowledge the anonymous referee for the detailed review and for the many
helpful suggestions which allowed us to improve the manuscript. This work was
supported by the Brazilian CNPq, FAPESP, and CAPES. This work has made use of
the facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul),
funded by FAPESP (grant 2009/54006-4) and INCT-A. This work has made use of
data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Amôres et al. (2017) Amôres, E. B., Robin, A. C., & Reylé, C. 2017, A&A, 602, A67
  * Antoja et al. (2018) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360
  * Barros et al. (2016) Barros, D. A., Lépine, J. R. D., & Dias, W. S. 2016, A&A, 593, A108
  * Barros et al. (2020) Barros, D. A., Pérez-Villegas, A., Lépine, J. R. D., Michtchenko, T. A., & Vieira, R. S. S. 2020, ApJ, 888, 75
  * Barros et al. (2021) Barros, D. A., Pérez-Villegas, A., Michtchenko, T. A., & Lépine, J. R. D. 2021, Frontiers in Astronomy and Space Sciences, 8, 48
  * Binney & Schönrich (2018) Binney, J. & Schönrich, R. 2018, MNRAS, 481, 1501
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bland-Hawthorn et al. (2019) Bland-Hawthorn, J., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 486, 1167
  * Bovy & Rix (2013) Bovy, J. & Rix, H.-W. 2013, ApJ, 779, 115
  * Burton & Gordon (1978) Burton, W. B. & Gordon, M. A. 1978, A&A, 63, 7
  * Carlin et al. (2013) Carlin, J. L., DeLaunay, J., Newberg, H. J., et al. 2013, ApJ, 777, L5
  * Clemens (1985) Clemens, D. P. 1985, ApJ, 295, 422
  * Darling & Widrow (2019) Darling, K. & Widrow, L. M. 2019, MNRAS, 484, 1050
  * Drimmel & Poggio (2018) Drimmel, R. & Poggio, E. 2018, Research Notes of the American Astronomical Society, 2, 210
  * Faure et al. (2014) Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564
  * Feldmann & Spolyar (2015) Feldmann, R. & Spolyar, D. 2015, MNRAS, 446, 1000
  * Fich et al. (1989) Fich, M., Blitz, L., & Stark, A. A. 1989, ApJ, 342, 272
  * Fragkoudi et al. (2019) Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 616, A11
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Ghosh et al. (2022) Ghosh, S., Debattista, V. P., & Khachaturyants, T. 2022, MNRAS, 511, 784
  * Gómez et al. (2013) Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  * Hunt et al. (2018) Hunt, J. A. S., Hong, J., Bovy, J., Kawata, D., & Grand, R. J. J. 2018, MNRAS, 481, 3794
  * Hunt et al. (2022) Hunt, J. A. S., Price-Whelan, A. M., Johnston, K. V., & Darragh-Ford, E. 2022, MNRAS, 516, L7
  * Kawata et al. (2018) Kawata, D., Baba, J., Ciucǎ, I., et al. 2018, MNRAS, 479, L108
  * Khachaturyants et al. (2022) Khachaturyants, T., Debattista, V. P., Ghosh, S., Beraldo e Silva, L., & Daniel, K. J. 2022, MNRAS, 517, L55
  * Khanna et al. (2019) Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 489, 4962
  * Khoperskov et al. (2019) Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6
  * Laporte et al. (2019) Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134
  * López-Corredoira et al. (2002) López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., & Hammersley, P. L. 2002, A&A, 394, 883
  * Michtchenko et al. (2019) Michtchenko, T. A., Barros, D. A., Pérez-Villegas, A., & Lépine, J. R. D. 2019, ApJ, 876, 36
  * Michtchenko et al. (2002) Michtchenko, T. A., Lazzaro, D., Ferraz-Mello, S., & Roig, F. 2002, Icarus, 158, 343
  * Michtchenko et al. (2018) Michtchenko, T. A., Lépine, J. R. D., Pérez-Villegas, A., Vieira, R. S. S., & Barros, D. A. 2018, ApJ, 863, L37
  * Michtchenko et al. (2017) Michtchenko, T. A., Vieira, R. S. S., Barros, D. A., & Lépine, J. R. D. 2017, A&A, 597, A39
  * Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  * Momany et al. (2006) Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515
  * Monari et al. (2015) Monari, G., Famaey, B., & Siebert, A. 2015, MNRAS, 452, 747
  * Rastorguev et al. (2017) Rastorguev, A. S., Utkin, N. D., Zabolotskikh, M. V., et al. 2017, Astrophysical Bulletin, 72, 122
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
  * Robin et al. (2022) Robin, A. C., Bienaymé, O., Salomon, J. B., et al. 2022, A&A, 667, A98
  * Smith et al. (2015) Smith, R., Flynn, C., Candlish, G. N., Fellhauer, M., & Gibson, B. K. 2015, MNRAS, 448, 2934
  * Wang et al. (2020) Wang, H. F., Huang, Y., Zhang, H. W., et al. 2020, ApJ, 902, 70
  * Wang et al. (2018) Wang, H.-F., Liu, C., Xu, Y., Wan, J.-C., & Deng, L. 2018, MNRAS, 478, 3367
  * Widrow et al. (2012) Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41
  * Williams et al. (2013) Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101


# Immunocto: a massive immune cell database auto-generated for histopathology

Mikaël Simard
Medical Physics, UCL, London, UK
<EMAIL_ADDRESS>
Zhuoyan Shen
Medical Physics, UCL, London, UK
Maria A. Hawkins
Medical Physics, UCL, London, UK
Radiotherapy, UCLH, London, UK
Charles-Antoine Collins-Fekete
Medical Physics, UCL, London, UK

###### Abstract

With the advent of novel cancer treatment options such as immunotherapy,
studying the tumour immune micro-environment is crucial to inform on prognosis
and understand response to therapeutic agents. A key approach to
characterising the tumour immune micro-environment may be through combining
(1) digitised microscopic high-resolution optical images of hematoxylin and
eosin (H&E) stained tissue sections obtained in routine histopathology
examinations with (2) automated immune cell detection and classification
methods. However, current individual immune cell classification models for
digital pathology present relatively poor performance. This is mainly due to
the limited size of currently available datasets of individual immune cells, a
consequence of the time-consuming and difficult problem of manually annotating
immune cells on digitised H&E whole slide images. In that context, we
introduce Immunocto, a massive, multi-million automatically generated database
of 6,848,454 human cells, including 2,282,818 immune cells distributed across
4 subtypes: CD4+ T cell lymphocytes, CD8+ T cell lymphocytes, B cell
lymphocytes, and macrophages. For each cell, we provide a $64\times 64$
pixels2 H&E image at $\mathbf{40}\times$ magnification, along with a binary
mask of the nucleus and a label. To create Immunocto, we combined open-source
models and data to automatically generate the majority of contours and labels.
The cells are obtained from a matched H&E and immunofluorescence colorectal
dataset from the Orion platform, while contours are obtained using the Segment
Anything Model. An advantage of Immunocto is that immunofluorescence data
largely facilitates cell type recognition and thus requires limited
annotations. A classifier trained on H&E images from Immunocto produces an
average $F_{1}$ score of 0.74 to differentiate the 4 immune cell subtypes and
other cells. For the more general task of identifying lymphocytes, the
classifier trained on Immunocto strongly outperforms a state of the art cell
segmentation and classification model, producing an $F_{1}$ score of 0.84
compared to 0.10. Immunocto can be downloaded at:
https://zenodo.org/uploads/11073373.

##  1 Introduction

###  1.1 Overview

Histopathology relates to the observation of abnormal tissue at the
microscopic scale. The gold standard methodology for histopathology is the
study of fixed tissue sections under a microscope stained with haematoxylin
and eosin (H&E). To identify diseases, pattern recognition at various scales
on tissue samples is performed by a pathologist. The advent of digital
pathology, where tissue sections are digitised to create gigapixel whole slide
images (WSIs), opens the possibility of creating large databases to train
models that may aid pathologists in their study of diseases and decision
making processes.

In the study of cancer, understanding the role of the tumour immune micro-
environment (TIME) is essential [1], as it can advise on prognosis, guide
treatment [2, 3], or be examined retrospectively to understand response to
therapies. For instance, an exploration of the TIME may reveal biomarkers that
could be useful in identifying patient populations that respond well to
immunotherapy [1].

Microscopic high-resolution imaging, where immune cell location and function
can be identified, is a critical next step to understand the TIME, as it can
provide information on cellular proportions, heterogeneity and spatial
distribution [1]. However, manual annotation of immune cells by pathologists
is time consuming and identifying various immune cell subtypes can be tedious
on H&E WSIs [4]. As such, there is growing interest in developing models that
can automatically segment and classify immune cells on WSIs [5].

As a first step towards building models that can recognise immune cells in
WSIs, we introduce Immunocto, a high-resolution ($\mathbf{40}\times$
magnification, pixel size of 0.325 $\mu$m) massive database of 2,282,818
immune cells distributed across 4 immune cell subtypes (CD4+ T cells, CD8+ T
cells, B cells, and macrophages). To our knowledge, Immunocto is the largest
available dataset of immune cells extracted from H&E WSIs by an order of
magnitude (table 1).

The main advantage of the methodology used to create the database is that it
requires minimal manual annotations, as it combines the Segment Anything Model
(SAM) [6] to extract candidate cells and relies on fully registered
immunofluorescence (IF) data from the Orion platform [7] to provide labels for
candidate cells.

Table 1: Fully labelled and segmented immune cells datasets ($\geq$ 10000 cells) on H&E WSIs. For cell subtypes, L is lymphocyte, M is macrophage, N is neutrophil, P is plasma, and E is eosinophil. Dataset | Cells | Cancer location | Cell subtypes
---|---|---|---
MoNuSAC [8] | 25,157 | Breast, Kidney, Lung, Prostate | L, M, N
PanNuke [9] | 32,276 | 19 sites | L, M, P
NuCLS [10] | 47,113 | Breast | E, L, M, N, P
Lizard [4] | 138,307 | 19 sites | E, L, N, P
Immunocto | 2,282,818 | Colon | L (CD4+ T cell, CD8+ T cell, B cell), M

###  1.2 Existing work

There have been considerable efforts in building labelled datasets for nuclear
segmentation and/or classification on H&E WSIs. The MoNoSAC dataset [8]
contains over 25,000 manually labelled immune cells. Manual labelling of
various type of immune cells, which often present similar morphological
features, can be cumbersome, requires extensive pathologist time and does not
scale well. The PanNuke database [9] rather relied on active learning, where
pre-labelled data is used to train a cell subtype classifier, and the results
are manually verified by pathologists. The PanNuke database is obtained after
7 iterations between model training and pathologist verification. Still,
manual investigation has a limited inter-observer reproducibility due to the
above-mentioned reason, which limits the accuracy of this methodology.
Alternatively, the larger NuCLS database [10] was obtained with a
crowdsourcing approach combining pathologists and non-pathologists (medical
students and graduates), which also needs important human intervention and is
still limited in accuracy by human identification on H&E images only. More
recently, the Lizard database [4] uses a four-step approach - iterative
segmentation and classification, cell boundary refinement, inflammatory nuclei
subtyping, and manual class refinement. Each step involves important
pathologist annotation, and the overall process includes the manual annotation
of over 70,000 immune cells.

There are two main limitations of the above datasets which motivated the
development of Immunocto. First, existing datasets involve important human
intervention to generate cell labels and/or nuclei masks, which limits the
number of labelled immune cells that can be produced. Second, as the above
datasets are obtained using H&E data only, the accuracy of the underlying
labels is limited. For instance, the Lizard database [4] has reported limited
inter-observer agreement for immune cell identification.

Furthermore, there have only been few investigations into subtyping
lymphocytes, which is a difficult and subjective task on H&E images due to
morphological similarities between lymphocyte types [11, 12]. For instance,
lymphocytes can be, amongst others, classified as CD8+ killer (cytotoxic) T
cells, CD4+ helper T cells, or B cells. Each type has a different function,
and capturing the spatial distribution of lymphocytes is an important part of
understanding the TIME [1]. As an example, stratification of colorectal cancer
patients has been correlated with deep infiltration of CD8+ T cells [13], and
the CD4+/CD8+ ratio of tumour-infiltrating lymphocytes has high prognostic
value in triple-negative breast cancer [14].

###  1.3 Key Innovations

Immunocto presents two key innovations. First, we use the Segment Anything
Model (SAM) to generate high-quality object masks at the cellular level. This
includes, but is not limited to, nuclear structures, epithelial and
endothelial cells, muscular and fat cells, apoptotic bodies, blood cells,
immune cells, and a variety of gaps and holes in the tissue.

The second innovation lies in the usage of IF as a tool to facilitate cell
labelling. IF is a technique that uses antibodies conjugated with fluorophores
to detect specific proteins on or within cells. When incubated with tissue
sections, the fluorescent antibodies will selectively bind to specific
antigens. Irradiating the sample using a laser with a wavelength that matches
a fluorophore’s absorption spectrum will create an IF image of the tissue
specimen, revealing locations on the tissue where a particular antigen is
expressed. For instance, the presence of CD4+ lymphocytes can be revealed with
CD4 antibodies. In that context, IF can be used to characterise various cell
types based on their unique surface markers, including CD4+, CD8+ T cells, or
B cells, the latter with CD20 antibodies. When an IF image is co-localised
with an H$\&$E image, an opportunity arises to automatically create a large
scale single-cell database, considering the unique object identification
originating from the SAM.

In this study, we used the Orion dataset [7], (distributed under the MIT
license) which contains H&E images co-localised with 18 IF data channels, to
extract cells that are active on IF images associated with immune cells.

##  2 Creating Immunocto

The overall framework for creating Immunocto is illustrated in figure 1; the 5
steps shown in the figure are detailed in section 2.2. Briefly, SAM is used on
a single WSI to detect objects, which are mostly cell nuclei. Using co-
localised IF data and a thresholding approach, 20,000 candidate immune cells
are automatically extracted from the SAM-derived masks. The cells are manually
reviewed by non-pathologists by looking at relevant IF channels. The result is
Immunocto $V_{0}$, the first iteration of the database containing 11,321
immune cells and 9048 other cells. This database is used to train a model
which predicts an immune cell type using H&E, IF images, and the binary mask
from SAM. The model is ran on 40 WSIs with matched H&E and IF data, and a
final, highly specific database, Immunocto $V_{1}$, is extracted by only
preserving cells with high classification probability (> 95%).

Figure 1: Workflow to create the Immunocto database.

###  2.1 Raw data source

To create Immunocto, we use 40 colorectal cancer patients with publically
available, registered H&E and IF data from the Orion platform proposed by Lin
et al. [7]. WSIs are obtained at $\mathbf{40}\times$ magnification with a
pixel size of 0.325 $\mu$m. IF data refers to a set of 18 fluorescence images
which are registered to H&E whole slide images and obtained via incubation
with various IF markers rather than H&E staining. Briefly, a tissue section is
incubated simultaneously with 18 IF markers, each tagging a specific type of
cell or structure. For instance, the Hoechst marker identifies nuclei, while
the CD20 marker selectively tags B cell lymphocytes. A fluorescence microscopy
system using various excitation lasers and tunable filters captures an image
for each IF marker, and the pixel-wise intensity of each image is proportional
to the concentration of the dye in that pixel. H&E staining and imaging is
then performed on the same tissue section.

Examples of H&E along with IF data can be seen in figure 1, as well as figure
3 of the supplementary material. With simple thresholding on selected IF
channels, one can extract regions of interest associated with specific cell
types and generate ground truth labels.

###  2.2 Database creation steps

#### Detecting cells with Segment Anything.

Cell nuclei in a given WSI can be accurately segmented using SAM [6] (ViT-H
architecture) without further fine tuning, provided adequate subsequent mask
filtering. To extract nuclei masks, a WSI is tessellated into non-overlapping
tiles of 512$\times$512 pixels2. Each tile is resampled to 1024$\times$1024
pixels2, and a grid of 64$\times$64 uniformly distributed points across the
image is used as a point prompt. Resulting masks are filtered with the
following empirically defined criteria: predicted intersection over union
(IoU) score of > 0.81, a stability score > 0.85, and an area between 42 and
4320 pixels2 to cover typical immune cell sizes. Masks with overlapping
bounding boxes of IoU larger than 0.1 are discarded via non-maximum
suppression using the `batched_nms` operator from PyTorch, with scores
interpreted as the original IoU provided by SAM. Resulting masks are arrays of
64$\times$64 pixels2. An example of cell detection with SAM over a large
region of interest is shown in the first three panels of figure 2.

#### Generating immune cell candidates via automatic thresholding.

Applying SAM on a WSI results, on average, in >1,000,000 detected objects,
most of which are cell nuclei. Leveraging IF data, one can extract candidate
objects for immune cells using an intensity thresholding scheme detailed in
supplementary material section B. For instance, detected objects which present
high signal in the Hoechst, CD3e and CD4 IF channels are classified as CD4+ T
cell lymphocytes. Overall, 25,000 candidate cells are automatically generated
- 5000 for each of the four immune cell subtype, and 5000 other, non-immune
cells.

#### Manual reviewing using an expert panel.

The automatic thresholding method produces misclassified candidate cells due
to imperfections with the IF staining. This includes slight registration
errors between IF and H&E channels, non-specific targeting of IF markers and a
mismatch between the SAM masks, which contour the nuclei, and the IF signal,
typically located on the cytoplasm of cells around the nuclei. Therefore, all
candidate cells are manually reviewed by an expert panel (supplementary
material section B.2) and re-classified if needed. The presence of the IF
channels allows non-pathologists to classify cells, as classification is not
based on morphological features of H&E images but mainly on the intensity of
selected IF channels. Manual review was performed by three authors of this
study (MS, ZS, CACF) using Label Studio v1.10.1 (https://labelstud.io). Some
cells marked as uncertain by one individual were reviewed by at least another
person before assigning a label. The resulting database, Immunocto $V_{0}$,
contains 11,321 immune cells (2512 CD4+ T cells, 2999 CD8+ T cells, 2411 B
cells and 3399 macrophages), as well as 9048 other cells.

#### Training a immune cell classifier from H&E, IF and cell masks.

We hypothesise that the above steps - recognising immune cells based on
matched IF and H&E data - is a relatively simple task that can be automated
with a trained deep neural network. Using the initial dataset Immunocto
$V_{0}$, we trained a model that classifies a 12-channel 64$\times$64 pixels2
image (which combines H&E, selected IF channels and the binary mask from SAM)
into one of five cell types - CD4+ T cell, CD8+ T cell, B cell, macrophage,
and other, non-immune cell. All training details and classification
performance metrics are provided in supplementary material section C.

#### Automatic generation of a massive database of immune cells.

The model trained on joint H&E and IF data is applied to the entire cohort of
the Orion platform [7], which contains 40 colorectal cancer whole slide
images. There are overall 54 million objects detected with SAM on the 40 WSIs;
to maximise the quality of automatically generated labels, we only preserve
immune cells which have a probability (softmax of logit resulting from the
model) above 95%. This results in Immunocto $V_{1}$, the final Immunocto
dataset, which contains 2,282,818 immune cells (952,929 CD4+ T cells, 467,395
CD8+ T cells, 480,900 B cells and 381,594 macrophages). Additionally,
Immunocto $V_{1}$ includes 4,565,636 (twice the number of immune cells) other,
non-immune cells to provide a database that can be used to discriminate
generally between immune and non-immune cells. To ensure a wide range of
sampled cells, we randomly sampled the 4,565,636 other cells from the subset
of objects in all 40 WSIs that were classified as other cell with >90%
probability.

###  2.3 The final Immunocto database

The steps from section 2.2 result in Immunocto $V_{1}$, which is the final
Immunocto database containing a total of 6,848,454 cells. For each cell, we
provide (1) a $64\times 64$ pixels2 H&E image at $\mathbf{40}\times$
magnification, (2) a $64\times 64$ pixels2 binary mask of the nucleus, and (3)
a label. The central point of the bounding box surrounding the cell is located
coincides with the center of the images. We provide additional dataset
documentation in supplementary material section A.

##  3 Applications

###  3.1 Evaluation of H&E classifiers trained on Immunocto on an internal
hold-out dataset

To demonstrate the versatility and clinical relevance of the Immunocto
dataset, we trained a model that predicts the immune cell subtype (CD4+ T
cells, CD8+ T cells, B cells, macrophages or other, non-immune cell) using
solely H&E stained images and cell masks (which can be obtained from SAM, or
any other segmentation algorithm).

The approach aims to demonstrate the utility of Immunocto using data
accessible in routine clinical practice, where immunofluorescence is not
readily available. Such a model could prove invaluable in two key situations:
firstly, as a prior assessment of biopsies/resections to guide treatment
decisions. Secondly, as tool for analysing the TIME either in low-ressource
settings where advanced staining techniques may be inaccessible, or for
retrospective studies where only H&E data is available.

To evaluate the performance of the classifier based on H&E data, we set aside
an internal hold-out dataset consisting of a single patient from the Orion
dataset and used the remaining 39 patients to train the model; training
details are reported in supplementary material section D.

We assessed immune cell classification based on H&E images solely using two
different approaches. In the first approach, a classifier was trained using
patches of 64$\times$64 pixel2, which contains information about individual
cell morphology. When reporting results, this model is presented as morph.
(64$\times$64). In the second approach, the classifier was trained with larger
patches of 256$\times$256 pixel2, enabling it to integrate information
concerning the individual cell as well as its interactions with its immediate
surrounding environment - this model is presented as morph. + context
(256$\times$256). Results on the hold-out test dataset are shown in table 2.
The table showcases performance metrics including precision (P), recall (R),
$F_{1}$ score ($F_{1}$) for each cell type for classifiers with the two patch
sizes. We also reported aggregated scores for other tasks such as identifying
T cells (CD4+ and CD8+) and lymphocytes (T cells and B cells).

Table 2: Prediction metrics (precision P, recall R and $F_{1}$ score) on the Orion hold-out patient per cell type. Results are reported for each model - the morph. 64$\times$64 captures the cell morphology, while morph. + context 256$\times$256 additionally captures the interaction of the cell with its immediate environment. The last two cell types (T cells and lymphocytes) represent aggregated classes. | morph. (64$\times$64) | morph. + context (256$\times$256)
---|---|---
Cell Type | P/R | $F_{1}$ | P/R | $F_{1}$
B cells | 0.75/0.69 | 0.72 | 0.78/0.80 | 0.79
CD4+ T cells | 0.60/0.62 | 0.61 | 0.63/0.66 | 0.65
CD8+ T cells | 0.54/0.59 | 0.56 | 0.57/0.57 | 0.57
Macrophages | 0.61/0.89 | 0.73 | 0.65/0.94 | 0.77
Others | 0.99/0.93 | 0.96 | 0.99/0.94 | 0.96
T cells (CD4 $\&$ CD8) | 0.89/0.95 | 0.92 | 0.91/0.94 | 0.93
Lymph. (T $\&$ B cells) | 0.92/0.97 | 0.94 | 0.94/0.97 | 0.95

We find that classification of immune cells with the larger context window
maximises the $F_{1}$ score for all classification tasks shown in table 2. The
improvement in $F_{1}$ score is most important for immune cell subtyping
(notably B cells, CD4+ T cells and macrophages), while context only provides
limited improvements for general lymphocyte detection.

###  3.2 Comparison with a state-of-the-art cell detection model

The above classifiers are compared with the state-of-the-art cell detection
model HoVer-Net [5] trained on the extensive cell dataset Lizard [4], which
was the published database with the highest number of annotated immune cells
on H&E datasets prior to this work. HoVer-Net simultaneously performs cell
segmentation and classification on H&E images. Since images in the Lizard
database have a pixel size of 0.5$\mu$m and the hold-out Orion WSI has a pixel
size of 0.325$\mu$m, the hold-out WSI was downsampled to 0.5$\mu$m before
applying HoVer-Net.

As HoVer-Net was not trained for lymphocyte subtyping, all models are compared
on the binary task of identifying if a cell is a lymphocyte or not. For the
models trained on Immunocto, a cell is assumed to be a lymphocyte if it is
either a CD4+ T cell, CD8+ T cell or a B cell. The models were applied on two
large regions of interest (ROI) from the hold-out Orion patient. The first
region, ROI 1, is a 3849$\times$3850 pixels2 patch containing a small number
of lymphocytes, while the second, ROI 2, is a 2800$\times$3850 pixels2 patch
with a high lymphocyte count. Table 3 reports the quantitative performance of
the models on the two ROIs, while figure 2 illustrates an overlay of ROI 2
with the segmentation and classification results from the Immunocto models and
HoVer-Net. Figure 4 in supplementary material section E.2 illustrates the
overlay for ROI 1.

Table 3: Prediction metrics of lymphocyte detection on a region with sparse lymphocytes (ROI 1) and a region with dense lymphocytes (ROI 2) from the Orion hold-out patient.  ROI | # | # other | morph. (64$\times$64) | morph. + context (256$\times$256) | HoVer-Net
---|---|---|---|---|---
lymph. | cells | P/R | F1 | P/R | F1 | P/R | F1
1 | 635 | 10,369 | 0.43/0.83 | 0.56 | 0.45/0.86 | 0.59 | 0.53/0.20 | 0.29
2 | 4,724 | 4,566 | 0.85/0.93 | 0.89 | 0.81/0.96 | 0.88 | 0.54/0.03 | 0.06
Total | 5,359 | 14,935 | 0.77/0.92 | 0.84 | 0.75/0.96 | 0.84 | 0.54/0.05 | 0.10

Overall, we find that classifiers trained on Immunocto largely outperforms
HoVer-Net for lymphocyte identification. As noted in table 2, context does not
provide significant performance enhancements for the general task of
classifying lymphocytes. While the HoVer-Net model has generally appreciable
precision, its recall is limited. This can be observed in figure 2d, where
most of the lymphocytes could not be detected with HoVer-Net. The difference
in performance may be attributed to the large discrepancy in the number of
training examples for lymphocytes between the Immunocto and Lizard datasets.
In addition, the lymphocyte labels provided by Immunocto are based on IF and
may be more robust than the manual ones from Lizard, which have limited inter-
observer agreement, as noted in figure C.1 of [4].

Figure 2: Overlay of a H&E region of interest containing a high lymphocyte
count with predicted cell types from various models. (a) shows the ground
truth labels from the Immunocto database; (b) and (c) respectively show the
predictions from the morph. 64$\times$64 and morph. + context 256$\times$256
models trained on the remaining Immunocto data. (d) shows the predictions
using HoVer-Net.

###  3.3 Challenges and limitations

Lymphocyte subtyping using co-localised H&E and IF data is a relatively simple
computer vision classification task that results in high $F_{1}$ scores,
between 0.92-0.96 across classes (table 1 of supplementary material). However,
when classification is performed solely on H&E images, table 2 shows largely
reduced $F_{1}$ scores for lymphocyte differentiation. This suggests that
subtyping lymphocytes on H&E images only is a difficult task due to their
morphological similarities [11, 12]. Therefore, the Immunocto dataset could be
used to benchmark more advanced computer vision models for classification.

The main limitation of Immunocto is that it has been trained only on
colorectal cancer data, as opposed to most other databases reported in table
1; this may limit the generalisation capabilities of models trained on
Immunocto for pan-cancer immune cell detection.

##  4 Conclusion

Immunocto is the largest available database of segmented immune cells on H&E
whole slide images for histopathology. It is automatically generated using
matched immunofluorescence data and SAM for object segmentation. It contains
952,929 CD4+ T cells, 467,395 CD8+ T cells, 480,900 B cells and 381,594
macrophages, as well as 4,565,636 other non-immune cells. The database can be
used to train models to study the tumour immune micro-environment on routine
H&E whole slide images. We report state of the art performance on immune cell
detection and subtyping using simple classifiers trained on Immunocto.

Furthermore, the methodology used to generate Immunocto can be extended to
more cell types, given that matched IF data is available. The approach
demonstrates the benefits of using matched H&E and IF data to generate robust
databases for computational pathology.

## Acknowledgments and Disclosure of Funding

This project is supported by the UKRI Future Leaders Fellowship, No.
MR/T040785/1, the Radiation Research Unit at the Cancer Research UK City of
London Centre Award C7893/A28990, as well as the UKRI AI for Health Award
EP/Y020030/1. We acknowledge Rasha Abu-Eid for fruitful discussions.

## References

  * [1] Mikhail Binnewies, Edward W Roberts, Kelly Kersten, Vincent Chan, Douglas F Fearon, Miriam Merad, Lisa M Coussens, Dmitry I Gabrilovich, Suzanne Ostrand-Rosenberg, Catherine C Hedrick, et al.  Understanding the tumor immune microenvironment (time) for effective therapy.  Nature medicine, 24(5):541–550, 2018.
  * [2] Wolf Herman Fridman, Franck Pagès, Catherine Sautès-Fridman, and Jérôme Galon.  The immune contexture in human tumours: impact on clinical outcome.  Nature Reviews Cancer, 12(4):298–306, 2012.
  * [3] Valerie Chew, Han Chong Toh, and Jean-Pierre Abastado.  Immune microenvironment in tumor progression: characteristics and challenges for therapy.  Journal of oncology, 2012, 2012.
  * [4] Simon Graham, Mostafa Jahanifar, Ayesha Azam, Mohammed Nimir, Yee-Wah Tsang, Katherine Dodd, Emily Hero, Harvir Sahota, Atisha Tank, Ksenija Benes, et al.  Lizard: a large-scale dataset for colonic nuclear instance segmentation and classification.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 684–693, 2021.
  * [5] Simon Graham, Quoc Dang Vu, Shan E Ahmed Raza, Ayesha Azam, Yee Wah Tsang, Jin Tae Kwak, and Nasir Rajpoot.  Hover-net: Simultaneous segmentation and classification of nuclei in multi-tissue histology images.  Medical image analysis, 58:101563, 2019.
  * [6] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C. Berg, Wan-Yen Lo, Piotr Dollár, and Ross Girshick.  Segment anything, 2023.
  * [7] Jia-Ren Lin, Yu-An Chen, Daniel Campton, Jeremy Cooper, Shannon Coy, Clarence Yapp, Juliann B Tefft, Erin McCarty, Keith L Ligon, Scott J Rodig, et al.  High-plex immunofluorescence imaging and traditional histology of the same tissue section for discovering image-based biomarkers.  Nature cancer, 4(7):1036–1052, 2023.
  * [8] Ruchika Verma, Neeraj Kumar, Abhijeet Patil, Nikhil Cherian Kurian, Swapnil Rane, Simon Graham, Quoc Dang Vu, Mieke Zwager, Shan E Ahmed Raza, Nasir Rajpoot, et al.  Monusac2020: A multi-organ nuclei segmentation and classification challenge.  IEEE Transactions on Medical Imaging, 40(12):3413–3423, 2021.
  * [9] Jevgenij Gamper, Navid Alemi Koohbanani, Ksenija Benes, Simon Graham, Mostafa Jahanifar, Syed Ali Khurram, Ayesha Azam, Katherine Hewitt, and Nasir Rajpoot.  Pannuke dataset extension, insights and baselines.  arXiv preprint arXiv:2003.10778, 2020.
  * [10] Mohamed Amgad, Lamees A Atteya, Hagar Hussein, Kareem Hosny Mohammed, Ehab Hafiz, Maha AT Elsebaie, Ahmed M Alhusseiny, Mohamed Atef AlMoslemany, Abdelmagid M Elmatboly, Philip A Pappalardo, et al.  Nucls: A scalable crowdsourcing approach and dataset for nucleus classification and segmentation in breast cancer.  GigaScience, 11:giac037, 2022.
  * [11] Elaine L Alexander and Bruce Wetzel.  Human lymphocytes: similarity of b and t cell surface morphology.  Science, 188(4189):732–734, 1975.
  * [12] Dmitry I Strokotov, Maxim A Yurkin, Konstantin V Gilev, Dirk R Van Bockstaele, Alfons G Hoekstra, Nikolay B Rubtsov, and Valeri P Maltsev.  Is there a difference between t-and b-lymphocyte morphology?  Journal of biomedical optics, 14(6):064036–064036, 2009.
  * [13] Bernhard Mlecnik, Gabriela Bindea, Helen K Angell, Pauline Maby, Mihaela Angelova, David Tougeron, Sarah E Church, Lucie Lafontaine, Maria Fischer, Tessa Fredriksen, et al.  Integrative analyses of colorectal cancer show immunoscore is a stronger predictor of patient survival than microsatellite instability.  Immunity, 44(3):698–711, 2016.
  * [14] Kai Wang, Tiansheng Shen, Gene P Siegal, and Shi Wei.  The cd4/cd8 ratio of tumor-infiltrating lymphocytes at the tumor-host interface has prognostic value in triple-negative breast cancer.  Human pathology, 69:110–117, 2017.