# Reviving stochasticity: uncertainty in SMBH binary eccentricity is
unavoidable

Alexander Rawlings,1 Matias Mannerkoski,1 Peter H. Johansson,1 Thorsten Naab2
1 Department of Physics, Gustaf Hällströmin katu 2, FI-00014, University of
Helsinki, Finland
2 Max-Planck-Institut für Astrophysik, Karl-Schwarzchild-Str 1, D-85748
Garching, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study supermassive black hole (SMBH) binary eccentricity of equal-mass
galaxy mergers in $N$-body simulations with the Ketju code, which combines the
gadget-4 fast multipole gravity solver with accurate regularised integration
and post-Newtonian corrections around SMBHs. In simulations with realistic,
high eccentricity galactic merger orbits, the hard binary eccentricity is
found to be a non-linear function of the deflection angle in the SMBH orbit
during the final, nearly radial close encounter between the SMBHs before they
form a bound binary. This mapping between the deflection angle and the binary
eccentricity has no apparent resolution dependence in our simulations spanning
the resolution range of $1\times 10^{5}$–$8\times 10^{6}$ particles per
galaxy. The mapping is also captured using a simple model with an analytic
potential, indicating that it is driven by the interplay between a smooth
asymmetric stellar background potential and dynamical friction acting on the
SMBHs. Due to the non-linearity of this mapping, in eccentric major merger
configurations small, parsec-scale variations in the merger orbit can result
in binary eccentricities varying in nearly the full possible range between
$e=0$ and $e=1$. In idealised simulations, such variations are caused by
finite resolution effects, and convergence of the binary eccentricity can be
achieved with increasing resolution. However, in real galaxies, other
mechanisms such as nuclear gas and substructure that perturb the merger orbit
are likely to be significant enough for the binary eccentricity to be
effectively random. Our results indicate that the distribution of these
effectively random eccentricities can be studied using even moderate
resolution simulations.

###### keywords:

black hole physics – galaxies: kinematics and dynamics – methods: numerical –
software: simulations

††pubyear: 2023††pagerange: Reviving stochasticity: uncertainty in SMBH binary
eccentricity is unavoidable–Reviving stochasticity: uncertainty in SMBH binary
eccentricity is unavoidable

##  1 Introduction

Supermassive black holes (SMBHs) are believed to reside at the centres of all
massive galaxies (e.g. Kormendy & Ho, 2013). In the $\Lambda$CDM model as
galaxies grow through gas accretion and mergers (e.g. Volonteri et al., 2003;
Naab & Ostriker, 2017) their SMBHs are also expected to interact in a three-
phase merger process (Begelman et al., 1980).

Firstly, dynamical friction (DF, Chandrasekhar, 1943) acts to bring the SMBHs
from kiloparsec-scales down to parsec-scale separations, after which the SMBHs
form a bound binary with a semimajor axis $a$ and eccentricity $e$ (e.g.
Milosavljević & Merritt, 2001; Merritt, 2013). In the second phase, the SMBH
binary separation is reduced through sequential slingshot encounters with the
surrounding stellar distribution (Hills & Fullerton, 1980; Hills, 1983;
Quinlan, 1996; Rantala et al., 2018). Finally, at small subparsec separations
gravitational wave (GW) emission becomes the dominant mechanism by which the
SMBH binary can lose its remaining orbital energy and angular momentum, thus
driving the SMBHs to coalescence (Peters & Mathews, 1963; Peters, 1964).

The complex nature of SMBH coalescence in a galaxy merger setting necessitates
the use of numerical techniques (e.g. Berentzen et al., 2009; Khan et al.,
2011; Dosopoulou & Antonini, 2017; Mannerkoski et al., 2019) to provide
quantitative predictions for ongoing observational programmes such as ground-
based pulsar timing arrays (PTAs, Agazie et al., 2023; Antoniadis et al.,
2023; Xu et al., 2023; Zic et al., 2023) and the upcoming Laser Interferometer
Space Antenna (LISA, e.g. Amaro-Seoane et al., 2023). In particular, the
eccentricity of the binary significantly affects the GW emission, with higher
eccentricities resulting both in faster mergers as well as changes in the
emitted GW spectrum (Enoki & Nagashima, 2007; Huerta et al., 2015; Taylor et
al., 2016). Understanding SMBH binary eccentricity and how faithfully it is
captured in numerical simulations is thus critical for predicting and
interpreting observations done with instruments such as PTAs and LISA.

The SMBH binary merging process and its dependence on eccentricity and
resolution has been extensively studied using collisionless merger simulations
(e.g. Berentzen et al., 2009; Vasiliev et al., 2015; Bortolas et al., 2016;
Gualandris et al., 2017, 2022). Recently, Nasim et al. (2020) have argued that
the substantial scatter in SMBH binary eccentricity observed in gas-free
merger simulations is an artefact of poor phase space discretisation, and that
in the real Universe where SMBH masses are far greater than stellar masses
$(M_{\bullet}\gg m_{\star})$, SMBH binary eccentricity is a reliably
predictable quantity.

In this paper, we find that for realistic galaxy merger orbits the scatter in
the final SMBH binary eccentricity is due to the physical sensitivity of the
final eccentricity to small perturbations of the final nearly radial plunging
trajectory of the SMBHs before they become bound, even in the infinite
resolution limit. As a result, we argue that the stochasticity of the binary
eccentricity is an unavoidable physical feature of realistic galaxy mergers,
at least in the equal-mass case.

This paper is structured as follows. In Section 2 we briefly discuss the main
features of the Ketju code and our numerical simulations. In Section 3 we
study the eccentricity scatter in our simulations and also present a simple
model that captures the observed scatter. We discuss our results and their
implications in Section 4 and finally present our conclusions in Section 5.

Figure 1: Binary eccentricity $e$ (note the non-linear scale) and shifted time
$t^{\prime}=t-t_{\mathrm{bound}}$ for the $M_{\bullet}/m_{\star}=5000$
resolution simulations from the $e_{0}=0.90$ set (blue lines) and $e_{0}=0.99$
set (orange lines). The circle markers indicate the time when a hard binary
formed, with the preceding semitransparent lines indicating the transition
period when the hard binary is still forming, and the solid lines the hard
binary evolution. During the transition period, the Keplerian definition of
binary motion is affected by the intervening stellar mass, creating the
oscillatory artefacts in the semitransparent lines. In the $e_{0}=0.90$
simulation set, six simulations have a binary merger within
$50\,\mathrm{Myr}$.

##  2 Numerical Simulations

We construct a number of idealised galaxy merger simulations, which we evolve
with our new version of the Ketju code (Mannerkoski et al., 2023; Rantala et
al., 2017). The dynamics of SMBHs, and stars in a small region around them,
are integrated with an algorithmically regularised integrator (Rantala et al.,
2020), whereas the dynamics of the remaining particles is computed with the
gadget-4 (Springel et al., 2021) fast multiple method (FMM) with second order
multipoles. Together with hierarchical time integration this allows for
symmetric interactions and manifest momentum conservation. Ketju also includes
post-Newtonian (PN) correction terms up to order 3.5 between each pair of
SMBHs (Blanchet, 2014).

Our galaxy models represent the nuclear bulge of a gas-devoid elliptical
galaxy, and are designed to match exactly the models of Nasim et al. (2020).
Each galaxy consists of a stellar-only Dehnen (1993) sphere with shape
parameter $\gamma=0.5$ and scale radius $a=186\,\mathrm{pc}$, where the Dehnen
profile is given by:

$\rho_{\star}(r)=\frac{(3-\gamma)M_{\star}}{4\pi}\frac{a}{r^{\gamma}(r+a)^{(4-\gamma)}}.$
(1)

The total stellar mass is $M_{\star}=10^{10}\,M_{\sun}$, and at the centre of
the model galaxy a SMBH of mass $M_{\bullet}=10^{8}\,M_{\sun}$ is placed111We
set the length and mass scales of the system to match one of those presented
in Gualandris et al. (2022). However, as the performed simulations are
gravity-only, the system may be transformed to other mass and length scales
with the same $M_{\mathrm{tot}}/a$ ratio without affecting the results.. Our
galaxies lie on the observed $M_{\bullet}\text{--}\sigma$ relation presented
in van den Bosch (2016). We test seven different mass resolutions with a
varying number of stellar particles: $N_{\star}=\\{1.0\times 10^{5},2.5\times
10^{5},5.0\times 10^{5},1.0\times 10^{6},2.0\times 10^{6},4.0\times
10^{6},8.0\times 10^{6}\\}$, corresponding to
$M_{\bullet}/m_{\star}=1000\textnormal{--}80000$.

We then construct isolated merger initial conditions by placing two galaxies
on two different elliptical orbits with eccentricities of $e_{0}=0.90$ and
$e_{0}=0.99$, at a fixed initial separation of $D=3.72\,\mathrm{kpc}$ and a
semimajor axis of $a_{0}=2.79\,\mathrm{kpc}$. The $e_{0}=0.90$ merger orbit
matches that used by Nasim et al. (2020), with the radial and tangential
velocities being consistent with the values reported in Gualandris et al.
(2022). For each orbital configuration, we run ten realisations for each mass
resolution, to account for stochasticity caused by the discretised phase
space. Interactions between stellar particles are softened with a softening
length of $\varepsilon=2.5\,\mathrm{pc}$, and the Ketju region radius is set
to $r_{\rm ketju}=3\varepsilon=7.5\,\mathrm{pc}$.

##  3 Results

Figure 2: Convergence of the eccentricity scatter $\sigma_{e}$ as a function
of mass resolution $M_{\bullet}/m_{\star}$ and the total number of stellar
particles $N_{\star,\mathrm{tot}}$. The results from Nasim et al. (2020) are
shown as teal triangles for comparison. The expected $1/\sqrt{N_{\star}}$
scaling is recovered for the $e_{0}=0.90$ mergers. For the $e_{0}=0.99$
mergers, $\sigma_{e}$ does not decrease until $M_{\bullet}/m_{\star}\geq
20000$.

###  3.1 Eccentricity scatter in simulations

The SMBH binary eccentricity for both the $e_{0}=0.90$ (blue) and $e_{0}=0.99$
(orange) orbits are shown in Figure 1 as a function of shifted time
$t^{\prime}=t-t_{\mathrm{bound}}$, where $t_{\mathrm{bound}}$ is the time when
the binary orbital energy $E$ becomes permanently negative. At the time
$t_{\mathrm{bound}}$ the SMBH binary is not yet isolated due to the large
quantity of intervening stellar mass, which causes the Keplerian definition of
eccentricity to oscillate and not instantaneously match exactly the actual
SMBH binary orbit. In Figure 1 we indicate the time of hard binary formation
by circle markers; from this point the Keplerian eccentricity appropriately
describes the binary orbit. The binary hardening radius $a_{\mathrm{h}}$ is
defined as (e.g. Merritt, 2006):

$a_{\mathrm{h}}=\frac{q}{(1+q)^{2}}\frac{r_{\mathrm{m}}}{4}=\frac{r_{\mathrm{m}}}{16},$
(2)

where $r_{\mathrm{m}}=r(m<2M_{\bullet,1})$ is the influence radius and $q$ is
the mass ratio between the SMBHs (for our simulations $q=1.0$). Before a hard
binary forms, both the $e_{0}=0.90$ and $e_{0}=0.99$ sets initially
demonstrate an overall decrease in eccentricity (semitransparent lines in
Figure 1) until $t^{\prime}\sim 1\,\mathrm{Myr}$. It should be noted that the
precise value of the eccentricity in this regime is not robust due to the
SMBHs orbiting in a potential still largely influenced by the stellar
background, however the overall trend is consistent with the binary orbital
geometry. After $t^{\prime}\sim 1\,\mathrm{Myr}$, some of the $e_{0}=0.99$
realisations continue to have a decreasing eccentricity, and all of the
$e_{0}=0.90$ realisations have an increasing eccentricity. We discuss the
mechanism for this in subsection 3.3.

The $e_{0}=0.99$ runs show a wide variation in eccentricity, where $e$ spans
almost the entire domain range $e=[0,1]$. As a test we also perform two sets
of ten runs using the $10^{6}$ particle resolution set up, but reduce the SMBH
mass to $5.0\times 10^{7}\,M_{\sun}$ and $1.0\times 10^{7}\,M_{\sun}$, and
observe the same qualitative spread in eccentricity as in the
$M_{\bullet}=10^{8}\,M_{\sun}$ case. Even though the $e_{0}=0.99$ runs have a
higher initial merger eccentricity than the $e_{0}=0.90$ runs, none of the
SMBH binaries in the shown $e_{0}=0.99$ set obtain high enough eccentricities
to undergo GW-induced coalescence during the $50\,\mathrm{Myr}$ timespan. The
$e_{0}=0.90$ set demonstrates six SMBH-binary mergers within
$50\,\mathrm{Myr}$ of forming a bound binary, seen as a rapid orbit
circularisation in Figure 1, which is captured self-consistently using Ketju.

To characterise the scatter in eccentricity, we determine the mean
eccentricity $e_{\mathrm{h}}$ over five orbital periods centred on the orbit
within which the SMBH binary has become hard (equation (2)). We characterise
the inter-simulation eccentricity scatter in the mean values of
$e_{\mathrm{h}}$ with the standard deviation, denoted as $\sigma_{e}$.

We show the dependence of $\sigma_{e}$ on mass resolution for the $e_{0}=0.90$
orbit in Figure 2 with blue circle markers, and for the $e_{0}=0.99$ orbit
using orange square markers.

For the $e_{0}=0.90$ orbit, the convergence of $\sigma_{e}$ scales as
$1/\sqrt{N_{\star,\mathrm{tot}}}$, where $N_{\star,\mathrm{tot}}$ is the total
number of stellar particles in the merger. The values of $\sigma_{e}$ we
report are quantitatively similar to the values found by Nasim et al. (2020)
for the same system, as can be seen by the teal triangles in Figure 2. The
variation in eccentricity does not show the same scaling for the $e_{0}=0.99$
orbit as the $e_{0}=0.90$ orbit. For mass resolutions
$M_{\bullet}/m_{\star}\leq 10000$, the value of $\sigma_{e}$ is almost
constant, before significantly dropping at higher mass resolutions.

Figure 3: Orbits from $t=0\,\mathrm{Myr}$ to a time shortly after
$t_{\mathrm{bound}}$ of a single SMBH in two representative realisations of
the $e_{0}=0.90$ mergers (left) and the $e_{0}=0.99$ mergers (right), which by
symmetry of the equal mass system is a reflection of the second BH orbit about
$x=-z$. The line gradient shows the shifted time
$t^{\prime}=t-t_{\mathrm{bound}}$, however the colouring for
$|t^{\prime}|>2\,\mathrm{Myr}$ is constant for visual clarity. In each inset
panel, $\theta_{\mathrm{defl}}$ is the deflection angle at the pericentre
between the arrows. The circle marker indicates the SMBH position when a bound
binary is formed.

###  3.2 Binary binding process and the scatter in eccentricity

To understand the observed scatter in eccentricity, we investigate the binary
binding process in each simulation. Before the SMBHs are bound, the primary
influence on the motion is from the galactic merger potential; two
realisations for both merger orbits are shown in Figure 3. The initial stages
of the orbit show indiscernible variation between realisations, however after
a particularly strong interaction between the SMBHs during a pericentre
passage, the orbits deviate between realisations (inset panels, Figure 3), and
evolve to different binary eccentricities.

To quantify the strength of the hard-scattering process that randomises the
SMBH orbit, we measure the effective two-body deflection angle, defined as:

$\theta_{\mathrm{defl}}=2\arctan\left(\frac{GM}{L\sqrt{2E}}\right)=2\arctan\left(\frac{b_{90}}{b}\right)$
(3)

where $M=M_{\bullet,1}+M_{\bullet,2}$, $L$ is the magnitude of the SMBH system
angular momentum vector, $E$ the orbital energy at the time of the pericentre
passage, $b$ the effective impact parameter, and $b_{90}$ the $90\degr$
deflection radius (Binney & Tremaine, 2008).

We observe clear evidence for a relationship between the deflection angle
$\theta_{\mathrm{defl}}$ and the resulting eccentricity at the time the SMBH
binary becomes hard, as shown with the data points in Figure 4. We also
observe a similar relationship in our two test simulation sets with the
reduced SMBH masses.

Figure 4:  Deflection angle $\theta_{\mathrm{defl}}$ and resulting hard binary
eccentricity in the $e_{0}=0.90$ (left) and $e_{0}=0.99$ (right) simulations,
compared to fitted model curves. The solid lines show the results from the
analytic model with $e_{\mathrm{s}}=0.91$, while the semitransparent lines
show the results for $e_{\mathrm{s}}=0.90$ and $e_{\mathrm{s}}=0.92$. The
$e_{0}=0.90$ model curves have been shifted left by $14\degr$ to better match
the simulation data, with the shift indicated by an arrow. The marginal
histograms show kernel density estimates of the simulation data. The inset
panel for the $e_{0}=0.90$ data shows a zoom-in with
$\theta_{\mathrm{defl}}=[52\degr,73\degr]$ and $e_{\mathrm{h}}=[0.96,1.00]$.

###  3.3 Reproducing the eccentricity behaviour with a simple model

The dependency of the binary eccentricity on the deflection angle
$\theta_{\mathrm{defl}}$ during the binary formation can be reproduced using a
simple analytic model, which includes only the essential components relevant
to this process. Taking the orbit to lie in the $xy$-plane, the relative
motion of the two equal-mass SMBHs in this model follows the equation of
motion

$\ddot{\bm{x}}=-\frac{2GM_{\bullet}}{|\bm{x}|^{3}}\bm{x}+\bm{a}_{\mathrm{bg}}+\bm{a}_{\mathrm{DF}},$
(4)

where $\bm{x}=(x,y)$ is the separation vector of the SMBHs,
$M_{\bullet}=10^{8}\,M_{\sun}$ is the mass of a single SMBH,
$\bm{a}_{\mathrm{bg}}$ the acceleration due to the asymmetric background
potential, and $\bm{a}_{\mathrm{DF}}$ is the acceleration due to DF.

The stellar background is modelled as a constant density spheroidal potential
(e.g. Binney & Tremaine, 2008):

$\Phi_{\mathrm{bg}}(\bm{x})=\pi G\rho(A_{x}x^{2}+A_{y}y^{2}),$ (5)

with the acceleration given by

$\bm{a}_{\mathrm{bg}}(\bm{x})=-\nabla\Phi_{\mathrm{bg}}(\bm{x}).$ (6)

The $A$ coefficients are related to the eccentricity $e_{\mathrm{s}}$ of the
spheroid as

$\displaystyle A_{x}$
$\displaystyle=2\left(\frac{1-e_{\mathrm{s}}^{2}}{e_{\mathrm{s}}^{2}}\right)\left[\frac{1}{2e_{\mathrm{s}}}\ln\left(\frac{1+e_{\mathrm{s}}}{1-e_{\mathrm{s}}}\right)-1\right]$
(7) $\displaystyle A_{y}$
$\displaystyle=\frac{1-e_{\mathrm{s}}^{2}}{e_{\mathrm{s}}^{2}}\left[\frac{1}{1-e_{\mathrm{s}}^{2}}-\frac{1}{2e_{\mathrm{s}}}\ln\left(\frac{1+e_{\mathrm{s}}}{1-e_{\mathrm{s}}}\right)\right],$
(8)

with the long axis aligned along the $x$-axis. The stellar density is set to
$\rho=300\,M_{\sun}\,\mathrm{pc}^{-3}$, which approximately matches the values
seen in the simulations when the SMBH binary is becoming bound.

Figure 5:  Top: Sample orbits computed with the analytic model for different
background potential eccentricities $e_{\mathrm{s}}$, with otherwise identical
initial conditions. The dashed lines show the isopotential contours of the
model. Only the early part of the orbital evolution is shown for clarity.
Bottom: Orbital eccentricity $e$ of the model orbits (note the non-linear
scale). The vertical line marks the end of the period shown in the top panels.

The DF acting on a single BH is modelled using the Chandrasekhar (1943)
formula assuming a Maxwellian distribution with a constant velocity dispersion
$\sigma_{\star}$ (e.g. Binney & Tremaine, 2008):

$\bm{a}_{1,\mathrm{DF}}=-\frac{4\pi
G^{2}M_{\bullet}\rho\ln{\Lambda}}{|\bm{v}_{1}|^{3}}\left(\operatorname{erf}(X)-\frac{2X}{\sqrt{\pi}}\exp(-X^{2})\right)\bm{v}_{1},$
(9)

where $\bm{v}_{1}=\dot{\bm{x}}/2$ is the velocity of a single BH, and
$X=|\bm{v}_{1}|/(\sqrt{2}\sigma_{\star})$. We set
${\sigma_{\star}=200\,\mathrm{km\,s^{-1}}}$ based on the value measured from
the simulations. The value of the Coulomb logarithm is expected to be in the
range $\ln{\Lambda}\sim 4\textnormal{--}5$ based on the size of the stellar
system, and we find that $\ln{\Lambda}=4.7$ gives results that agree well with
the simulations.

The effect of DF weakens as the SMBH binary orbit shrinks. To account for
this, when the binary is bound we multiply the acceleration given by equation
(9) with the smooth cut-off function

$f(a)=\frac{1}{1+\exp[(a_{\mathrm{c}}-a)/d_{\mathrm{c}}]}.$ (10)

Here $a$ is the semimajor axis of a bound binary, and the cut-off scales are
$a_{\mathrm{c}}=2a_{\mathrm{h}}$ and $d_{\mathrm{c}}=0.5a_{\mathrm{h}}$, with
$a_{\mathrm{h}}$ defined by equation (2). The total DF term is then

$\bm{a}_{\mathrm{DF}}=2f(a)\bm{a}_{1,\mathrm{DF}}.$ (11)

To mimic the nearly linearly plunging orbits before the scattering event seen
in Figure 3, we specify the initial conditions as

$\displaystyle\bm{x}_{0}$ $\displaystyle=(25\,\mathrm{pc},b)$ (12)
$\displaystyle\dot{\bm{x}}_{0}$ $\displaystyle=(-v_{0},0).$ (13)

In order to match the typical SMBH relative velocity seen at this separation
in the simulations we set $v_{0}=450\,\mathrm{km\,s^{-1}}$ for the
$e_{0}=0.90$ case, and $v_{0}=560\,\mathrm{km\,s^{-1}}$ for the the
$e_{0}=0.99$ case.

The spheroid eccentricity parameter $e_{\mathrm{s}}$ is not easily measured
from the simulations, due to the presence of stellar components that remain
tightly bound to the SMBHs and obscure the background potential relevant for
the dynamics, in combination with a constantly evolving potential. However,
only sufficiently large values of $e_{\mathrm{s}}$ allow for low binary
eccentricities to be produced, as is shown by Figure 5. Performing
calculations with different values of $e_{\mathrm{s}}$, we found a good fit to
the simulation results for $e_{\mathrm{s}}\approx 0.9$.

We compute the resulting eccentricities for a range of impact parameters up to
$b=20\,\mathrm{pc}$ for the two initial velocities, by solving the equation of
motion (4) until the binary has become hard using the error controlled 8th
order Runge–Kutta method DOP853 included in the SciPy library (Virtanen et
al., 2020). For reference, the $90\degr$ deflection radii are $b_{90}\approx
5.8\,\mathrm{pc}$ and $b_{90}\approx 3.3\,\mathrm{pc}$ for the $e_{0}=0.90$
and $e_{0}=0.99$ cases, respectively.

The resulting model curves are shown together with the simulation data in
Figure 4. To correctly match the data in the $e_{0}=0.90$ case, the model
curve has been shifted to the left by $14\degr$. This shift is required likely
due to the rotation of the stellar component tilting the background potential
relative to the SMBH trajectories in the simulations, with the effect being
only evident in the $e_{0}=0.90$ case due to the less radial merger orbit. In
the $e_{0}=0.99$ case we can also see that the simulation data has significant
scatter around the model curve in the $\theta_{\mathrm{defl}}\gtrsim 120\degr$
region, although the simple model does appear to capture the mean behaviour
relatively well, apart from the very largest values of
$\theta_{\mathrm{defl}}$. The behaviour of the model in this region is also
sensitive to the exact values of $e_{\mathrm{s}}$ and $\ln{\Lambda}$, which
might also be related to the large scatter seen in the simulation data.

However, in general the behaviour of the simulation data is captured well,
with the model correctly producing the two main eccentricity minima, as well
as the $e\approx 1$ region between them. The predicted minima occur for
trajectories where the torque from the background potential together with the
DF causes the SMBHs to loop around into a nearly circular orbit, resulting in
the low eccentricities seen in some of the $e_{0}=0.99$ simulation data.
Conversely, highly eccentric binaries are produced when the binary is trapped
into nearly radial oscillations along the main axes of the potential. These
mechanisms are contrasted in Figure 5.

##  4 Discussion

As shown by both the simulations and the analytic model, the eccentricity of
the hard SMBH binary in the studied merger configurations can span nearly the
full possible range in the interval $e_{\mathrm{h}}=[0,1]$, depending on
comparatively small, parsec-scale differences in the particular realisation of
the galactic merger orbit. On the other hand, in minor mergers where the
stellar background is less asymmetric when the SMBHs become bound, the scatter
in eccentricity is likely to be relatively low even in the case of a radial
merger orbit, since the system is less sensitive to the exact value of the
deflection angle, as seen in the $e_{\mathrm{s}}=0.2$ case in Figure 5.

Merger orbits with a lower initial eccentricity $e_{0}$ than the values in
this work can be expected to show less scatter in the binary eccentricity due
to the lack of a hard scattering event that is sensitive to slight
perturbations in the merger orbit, which was also seen in the study by
Gualandris et al. (2022). However, major mergers are expected to occur on
highly radial orbits based on cosmological simulations (e.g. Khochfar &
Burkert, 2006). In addition, the DF from the dark matter halo would in general
drive even lower eccentricity merger orbits to highly radial final plunges
when the initial separation of the galaxies is large enough.

The simulation setup used in this work is idealised, and covers only a small
part of the parameter space of merger configurations. However, similar
behaviour can be expected to occur also in more realistic models, since the
relevant dynamics occurs within a few hundred parsecs from the centre of the
merged galaxy, where we do not expect any significant effects from dark matter
or the outer parts of the galaxy, which are not included in our present
simulations. Preliminary results using steeper $(\gamma\gtrsim 1)$, more
observationally-motivated galaxy density profiles indicate that similar
behaviour occurs also in more realistic merger systems.

The variation of the deflection angle is caused by the random variation of the
impact parameter $b$ through equation (3). In numerical simulations with a
discretised phase space, the primary cause of random variation in $b$ is the
Brownian motion of the SMBHs before they form a bound binary (Merritt, 2001;
Bortolas et al., 2016). This effect scales as
$1/\sqrt{N_{\star,\mathrm{tot}}}$, which can explain the observed scaling of
$\sigma_{e}$ (Nasim et al., 2020) when the system falls into the region of
parameter space where the relation between $e$ and $b$ is approximately
linear. The deviation from this scaling seen in Figure 2 can then be explained
by the fact that the eccentricity is not globally a linear function of the
impact parameter. The non-linear mapping also explains why GW-induced mergers
were present within a short ($<50\,\mathrm{Myr}$) timeframe for the
$e_{0}=0.90$ mergers and not the $e_{0}=0.99$ mergers in Figure 1, and
highlights the difficulty in predicting the binary eccentricity
$e_{\mathrm{h}}$ from the initial merger orbit.

While the scatter in the impact parameter between different numerical
realisations of the same merger is due to relatively low number of particles
compared to real galaxies, we expect that uncertainty of a similar magnitude
is also present in real systems due to various mechanisms, such as
perturbations of the SMBH motion from gas and substructure in their stellar
environment, or larger-scale perturbations of the galactic merger orbit. In
addition, the SMBHs are also not necesssarily located exactly at the centre of
mass of the galaxy. For example, the work by Batcheldor et al. (2010) has
suggested that the SMBH in M87 may be displaced from the galactic photocentre
by up to $\sim 7\,\mathrm{pc}$ due to jet acceleration or gravitational recoil
kicks. Thus, it is not possible to give an exact prediction for the SMBH
binary eccentricity produced by a given merger, and instead we must focus on
predicting the distribution of eccentricities that can be produced by a given
merger configuration.

Predicting the eccentricity distribution requires the knowledge of the
distribution of impact parameters that can result from essentially identical
mergers as well as the mapping between the impact parameter, or deflection
angle, and the final eccentricity. Since the relation between the impact
parameter and the eccentricity does not appear to significantly depend on the
resolution, simulations at moderate mass resolutions of
$M_{\bullet}/m_{\star}\sim 10^{3}$ can be used to map out the relation by
performing a large number of merger simulations, possibly augmented by more
sophisticated versions of the analytic model presented here. The distribution
of impact parameters in real systems is however a much more difficult problem
to tackle, since in simulations numerical resolution effects are likely
dominant.

Even without such extensive studies on the full distribution of SMBH binary
eccentricities, the present results suggest that the distribution is likely to
be broad and contain a significant number of binaries at high eccentricities.
The possibility of a large fraction of high binary eccentricities has
implications in particular for studies of the expected gravitational wave
background (GWB) signal measured by PTAs, as in particular binaries with
$e_{\mathrm{h}}\gtrsim 0.9$ can retain significant eccentricities while their
GW emission is within the observable frequency band. This could be observable
in the shape of the GWB spectrum (Kelley et al., 2017), and would also affect
the likelihood of detecting individual SMBH binaries with PTAs (Kelley et al.,
2018). It is interesting to note that Bi et al. (2023) obtained a distribution
of eccentricities consistent with the expectation of a significant fraction of
SMBH binaries at high eccentricities based on analysis of the recent PTA
observations.

In the present work, the influences of rotation and SMBH orbital flips are not
considered. A counter-rotating stellar background with respect to the SMBH
binary angular momentum has a tendency to increase the binary eccentricity,
whereas a corotating stellar background tends to circularise the binary
(Sesana et al., 2011; Holley-Bockelmann & Khan, 2015). The SMBH angular
momentum vector can however undergo reversals as a result of gravitational
torques following pericentre passages of the binary (Rantala et al., 2019,
termed an orbital flip) an example of which is shown in the left inset of of
the left panel of Figure 3. The chaotic reversals of the binary angular
momentum act to randomise the relative rotation of the stellar background,
leading in general to complex eccentricity evolution (Nasim et al., 2021). A
result of a sustained co- or counter-rotating stellar background could be an
evolution of the binary eccentricity after $e_{\mathrm{h}}$ to a bimodal
distribution for a given merger configuration, if the stellar scattering
interactions were strong enough. Whilst these mechanisms are unlikely to alter
much our findings concerning the hard binary eccentricity, they may obfuscate
predictions for SMBH binary eccentricity prior to the merger, which is
critical for GW detection missions.

In contrast to the gas-free and fairly low density systems studied here, in
higher density or gas-rich systems the binary eccentricity can also evolve
significantly after the binary has become hard but before it has entered the
GW emission dominated regime. In particular, interactions with a circumbinary
accretion disc can lead to rapid evolution of the eccentricity if the
accretion rate is high. Studies of prograde coplanar circumbinary discs have
found that at almost all eccentricities the binaries evolve towards an
eccentricity of $e\approx 0.5$ (Zrake et al., 2021; D’Orazio & Duffell, 2021),
whereas retrograde discs lead to the eccentricity always increasing (Tiede &
D’Orazio, 2023). However, these studies are generally limited to
eccentricities below $e<0.9$. Our results indicate the need for extending such
studies to higher binary eccentricities, as well as to more general
configurations such as polar discs which become increasingly likely at high
eccentricities (e.g. Martin & Lubow, 2017).

##  5 Conclusions

In conclusion, we find that the variation in the hard binary eccentricity
$e_{\mathrm{h}}$ for a given system configuration is tightly correlated with
the deflection angle $\theta_{\mathrm{defl}}$, and thus the impact parameter
$b$. By using a simple, resolution-free analytic model of the SMBH scattering
process, we have demonstrated that uncertainty in SMBH binary eccentricity is
not caused solely by discretisation effects in galaxy merger simulations, but
is rather due to the physical sensitivity of the system to small changes in
the merger orbit, which can be caused by physical mechanisms in addition to
numerical discretisation effects. The results presented here justify extending
the investigation to more realistic galaxy merger scenarios, in order to
quantify the expected range of hard binary eccentricities, which is critical
for predictions for current and future gravitational wave observation
missions.

## acknowledgments

We thank the anonymous referee for their comments which contributed to
improving this work. We also thank Shihong Liao, Dimitrios Irodotou, and
Francesco Paolo Rizzuto for interesting discussions. A.R. acknowledges the
support by the University of Helsinki Research Foundation. A.R., M.M. and
P.H.J. acknowledge the support by the European Research Council via ERC
Consolidator Grant KETJU (no. 818930) and the support of the Academy of
Finland grant 339127. TN acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence
“ORIGINS”.

The numerical simulations used computational resources provided by the CSC –
IT centre for Science, Finland.

## Author contributions

We list here the roles and contributions of the authors according to the
Contributor Roles Taxonomy (CRediT). AR: Conceptualisation, Investigation,
Formal analysis, Data curation, Writing – original draft. MM:
Conceptualisation, Formal analysis, Writing – original draft. PHJ:
Supervision, Writing – review & editing. TN: Writing – review & editing.

## Software

Ketju (Mannerkoski et al., 2023; Rantala et al., 2017), gadget-4 (Springel et
al., 2021), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
Matplotlib (Hunter, 2007), pygad (Röttgers et al., 2020).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Agazie et al. (2023) Agazie G., et al., 2023, ApJ, 951, L8
  * Amaro-Seoane et al. (2023) Amaro-Seoane P., et al., 2023, Living Reviews in Relativity, 26, 2
  * Antoniadis et al. (2023) Antoniadis J., et al., 2023, arXiv e-prints, p. arXiv:2306.16227
  * Batcheldor et al. (2010) Batcheldor D., Robinson A., Axon D. J., Perlman E. S., Merritt D., 2010, ApJ, 717, L6
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Berentzen et al. (2009) Berentzen I., Preto M., Berczik P., Merritt D., Spurzem R., 2009, ApJ, 695, 455
  * Bi et al. (2023) Bi Y.-C., Wu Y.-M., Chen Z.-C., Huang Q.-G., 2023, arXiv e-prints, p. arXiv:2307.00722
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Blanchet (2014) Blanchet L., 2014, Living Reviews in Relativity, 17, 2
  * Bortolas et al. (2016) Bortolas E., Gualandris A., Dotti M., Spera M., Mapelli M., 2016, MNRAS, 461, 1023
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * D’Orazio & Duffell (2021) D’Orazio D. J., Duffell P. C., 2021, ApJ, 914, L21
  * Dehnen (1993) Dehnen W., 1993, MNRAS, 265, 250
  * Dosopoulou & Antonini (2017) Dosopoulou F., Antonini F., 2017, ApJ, 840, 31
  * Enoki & Nagashima (2007) Enoki M., Nagashima M., 2007, Progress of Theoretical Physics, 117, 241
  * Gualandris et al. (2017) Gualandris A., Read J. I., Dehnen W., Bortolas E., 2017, MNRAS, 464, 2301
  * Gualandris et al. (2022) Gualandris A., Khan F. M., Bortolas E., Bonetti M., Sesana A., Berczik P., Holley-Bockelmann K., 2022, MNRAS, 511, 4753
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hills (1983) Hills J. G., 1983, AJ, 88, 1269
  * Hills & Fullerton (1980) Hills J. G., Fullerton L. W., 1980, AJ, 85, 1281
  * Holley-Bockelmann & Khan (2015) Holley-Bockelmann K., Khan F. M., 2015, ApJ, 810, 139
  * Huerta et al. (2015) Huerta E. A., McWilliams S. T., Gair J. R., Taylor S. R., 2015, Phys. Rev. D, 92, 063010
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Kelley et al. (2017) Kelley L. Z., Blecha L., Hernquist L., Sesana A., Taylor S. R., 2017, MNRAS, 471, 4508
  * Kelley et al. (2018) Kelley L. Z., Blecha L., Hernquist L., Sesana A., Taylor S. R., 2018, MNRAS, 477, 964
  * Khan et al. (2011) Khan F. M., Just A., Merritt D., 2011, ApJ, 732, 89
  * Khochfar & Burkert (2006) Khochfar S., Burkert A., 2006, A&A, 445, 403
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Mannerkoski et al. (2019) Mannerkoski M., Johansson P. H., Pihajoki P., Rantala A., Naab T., 2019, ApJ, 887, 35
  * Mannerkoski et al. (2023) Mannerkoski M., Rawlings A., Johansson P. H., Naab T., Rantala A., Springel V., Irodotou D., Liao S., 2023, MNRAS,
  * Martin & Lubow (2017) Martin R. G., Lubow S. H., 2017, ApJ, 835, L28
  * Merritt (2001) Merritt D., 2001, ApJ, 556, 245
  * Merritt (2006) Merritt D., 2006, ApJ, 648, 976
  * Merritt (2013) Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei. Princeton University Press
  * Milosavljević & Merritt (2001) Milosavljević M., Merritt D., 2001, ApJ, 563, 34
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nasim et al. (2020) Nasim I., Gualandris A., Read J., Dehnen W., Delorme M., Antonini F., 2020, MNRAS, 497, 739
  * Nasim et al. (2021) Nasim I. T., Petrovich C., Nasim A., Dosopoulou F., Antonini F., 2021, MNRAS, 503, 498
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Quinlan (1996) Quinlan G. D., 1996, New Astron., 1, 35
  * Rantala et al. (2017) Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53
  * Rantala et al. (2018) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2018, ApJ, 864, 113
  * Rantala et al. (2019) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2019, ApJ, 872, L17
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Röttgers et al. (2020) Röttgers B., Naab T., Cernetic M., Davé R., Kauffmann G., Borthakur S., Foidl H., 2020, MNRAS, 496, 152
  * Sesana et al. (2011) Sesana A., Gualandris A., Dotti M., 2011, MNRAS, 415, L35
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Taylor et al. (2016) Taylor S. R., Huerta E. A., Gair J. R., McWilliams S. T., 2016, ApJ, 817, 70
  * Tiede & D’Orazio (2023) Tiede C., D’Orazio D. J., 2023, arXiv e-prints, p. arXiv:2307.03775
  * Vasiliev et al. (2015) Vasiliev E., Antonini F., Merritt D., 2015, ApJ, 810, 49
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Volonteri et al. (2003) Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
  * Xu et al. (2023) Xu H., et al., 2023, Research in Astronomy and Astrophysics, 23, 075024
  * Zic et al. (2023) Zic A., et al., 2023, arXiv e-prints, p. arXiv:2306.16230
  * Zrake et al. (2021) Zrake J., Tiede C., MacFadyen A., Haiman Z., 2021, ApJ, 909, L13
  * van den Bosch (2016) van den Bosch R. C. E., 2016, ApJ, 831, 134


# Dual Approaches to Express the Generalized Degree of Polarimetric Purity

A. Bhattacharya, S. Dey, A. C. Frery

##  1 Introduction

The degree of polarimetric purity is an invariant dimensionless quantity that
characterizes the closeness of a polarization state of a wave to a pure state
and is related to the Von Neumann entropy [1]. The polarimetric purity of a
plane wave characterized by the second-order statistics (i.e., the covariance
matrix) is uniquely described by the degree of polarization. However, the 2D
formalism is only applicable when the wave propagation direction is fixed.
This assumption is typical in optical and radar polarimetric measurements.
Therefore, one must consider all the components to describe the general state
of wave polarization. Starting from Samson [2], and Barakat [3], several
different concepts have been proposed in the literature to describe the 3D
degree of polarization [4, 5, 6, 7, 8, 9].

The generalized degree of polarimetric purity, $P_{nD}$ [10, 11, 12] for the
$n\times n$ Hermitian and positive semi-definite coherency matrix
$\mathbf{\Phi}$ is defined as,

$P_{nD}=\left\\{\frac{1}{n-1}\Bigg{[}\frac{n[\operatorname{tr}{(\mathbf{\Phi}^{2})}]}{(\operatorname{tr}{\mathbf{\Phi}})^{2}}-1\Bigg{]}\right\\}^{1/2}$
(1)

where $\operatorname{tr}(\mathbf{\Phi})$ is the trace of $\mathbf{\Phi}$. This
is an invariant dimensionless quantity satisfying, $0\leq P_{nD}\leq 1$. The
minimum value of the degree of purity $P_{nD}=0$ corresponds to a completely
random polarization state of a wave. In contrast, the maximum value $P_{nD}=1$
corresponds to a completely pure polarization state of a wave. In between the
two extremities, the wave is partially pure polarized.

In this work, we propose two distinct approaches to express the generalized
degree of polarimetric purity, $P_{nD}$. In the first approach, we utilize the
definition of the mean and standard deviation of real positive eigenvalues of
Hermitian positive semi-definite matrices [13]. In the second approach, we use
elementary concepts from vector calculus and align them with a matrix
decomposition procedure following certain notions from Lie algebra [14].
Finally, we establish the parity of the two approaches to compute the
generalized degree of polarimetric purity.

##  2 Approach I: Coefficient of Variation

The mean $(m)$ and the standard deviation $(s)$ of the real positive
eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}\geq 0$ for a
$n\times n$ coherency matrix $\mathbf{\Phi}$ are defined using a simple
function of the trace of the matrix and the trace of its square given in [15]
as,

$\displaystyle m$
$\displaystyle=\frac{1}{n}\operatorname{tr}(\mathbf{\Phi})=\frac{1}{n}\sum_{j=1}^{n}\lambda_{j},\text{
and}$ (2) $\displaystyle s^{2}$
$\displaystyle=\frac{1}{n}\Bigg{[}\sum_{j=1}^{n}\lambda_{j}^{2}-\frac{1}{n}\bigg{(}\sum_{j=1}^{n}{\lambda_{j}}\bigg{)}^{2}\Bigg{]}=\frac{\operatorname{tr}{(\mathbf{\Phi}^{2})}-(\operatorname{tr}{\mathbf{\Phi}})^{2}/n}{n},$
(3) $\displaystyle=\frac{\operatorname{tr}(\mathbf{\Phi}^{2})}{n}-m^{2}.$ (4)

Using these two quantities, we propose a new expression for the generalized
degree of polarimetric purity as,

$P_{nD}=\frac{s}{\sqrt{n-1}\,m}$ (5)

One can easily observe that for $n=2$, and $n=3$, the expressions for the
generalized degree of polarimetric purity can be related as [13],

$\displaystyle P_{2D}$
$\displaystyle=\left\\{\frac{2[\operatorname{tr}{(\mathbf{\Phi}^{2})}]}{(\operatorname{tr}{\mathbf{\Phi}})^{2}}-1\right\\}^{1/2}=s/m$
$\displaystyle P_{3D}$
$\displaystyle=\left\\{\frac{1}{2}\Bigg{[}\frac{3[\operatorname{tr}{(\mathbf{\Phi}^{2})}]}{(\operatorname{tr}{\mathbf{\Phi}})^{2}}-1\Bigg{]}\right\\}^{1/2}=s/\sqrt{2}m$

The proposed expression (5) is physically intuitive as it directly relates the
measure of polarimetric purity to the coefficient of variation (i.e., $s/m$)
of the eigenvalues of a $n\times n$ matrix.

Furthermore, let $\psi$ denote the angle between $\mathbf{\Phi}$ and the
identity matrix $\mathbf{I}_{n}$ in the space of $n\times n$ matrices.
Analogously, we can express this angle between the vector of eigenvalues
$\lambda_{1},\lambda_{2},\dots,\lambda_{n}$ and the equiangular line as [15],

$\psi=\cos^{-1}\left(\operatorname{tr}{\mathbf{\Phi}}/\sqrt{n\left(\operatorname{tr}{\mathbf{\Phi}}\right)^{2}}\right)$
(6)

therefore, using $\psi$ we can also express the generalized degree of
polarimetric purity as,

$P_{nD}=\frac{\tan{\psi}}{\sqrt{n-1}}$ (7)

which immediately relates, that $\tan{\psi}=s/m$.

##  3 Approach II: Direct Sum Decomposition

Let $\mathbf{\Phi}$ be decomposed as,

$\mathbf{\Phi}=\mathbf{\Phi}_{1}+\mathbf{\Phi}_{2}+\mathbf{\Phi}_{3}$ (8)

where,

$\displaystyle\mathbf{\Phi}_{1}$
$\displaystyle=\dfrac{(\mathbf{\Phi}+\mathbf{\Phi}^{*})}{2}-\dfrac{\operatorname{tr}(\mathbf{\Phi})}{n}\mathbf{I}_{n},\quad(\text{Traceless
symmetric matrix})$ (9) $\displaystyle\mathbf{\Phi}_{2}$
$\displaystyle=\dfrac{(\mathbf{\Phi}-\mathbf{\Phi}^{*})}{2},\quad(\text{Anti-
symmetric matrix})$ (10) $\displaystyle\mathbf{\Phi}_{3}$
$\displaystyle=\dfrac{\operatorname{tr}(\mathbf{\Phi})}{n}\mathbf{I}_{n},\quad(\text{Scalar
matrix})$ (11)

where $\mathbf{I}_{n}$ is the $n\times n$ identity matrix and
$\mathbf{\Phi}^{*}$ is the conjugate of $\mathbf{\Phi}$.

On the one hand, we can consider this categorization as a direct sum
decomposition of the Lie algebra $\mathfrak{gl}(n)$. It is known from the
literature [14] that the sub-algebra of traceless matrices is the Lie algebra
$\mathfrak{sl}(n)$ of the $SL(n)$ group (i.e., the special linear group). The
anti-symmetric matrices form the Lie algebra $\mathfrak{so}(n)$ of the $SO(n)$
group (i.e., the special orthogonal group).

On the other hand, we can interpret using elementary property from vector
calculus that the symmetric, trace-free derivative operation physically
relates to that of a shear [16]. Mathematically this operation is represented
by the matrix, $\mathbf{\Phi}_{1}$, which one can imagine as the gradient of a
vector field in an arbitrary direction (usually expressed by a Jacobian
matrix). However, the anti-symmetric matrix, $\mathbf{\Phi}_{2}$ represents
pure rotation (i.e., the curl operator).

At first, using the traceless symmetric matrix, $\mathbf{\Phi}_{1}$, we define
a quantity as,

$P_{ns}=\frac{\sqrt{\dfrac{n}{n-1}}{\left\lVert\mathbf{\Phi}_{1}\right\rVert}_{\text{F}}}{\operatorname{tr}(\mathbf{\Phi})}$
(12)

where ${\left\lVert\cdot\right\rVert}_{\text{F}}$ is the Frobenius norm of the
matrix. We call this quantity the degree of shearing. One can show that
$P_{ns}$ is invariant under unitary transformation. In particular, $P_{3s}$
can be considered as the degree of polarization of the real part of the
partially polarized $3\times 3$ intrinsic coherency matrix,
$Re(\mathbf{\Phi})$ [17]. Moreover, it is interesting to note that we can
relate the degree of shearing to the component of polarimetric purity (CPP)
proposed by Gil [18, 19] i.e., the degree of linear polarization, $P_{l}$, for
2D and 3D cases, and the degree of directionality, $P_{d}$, for the 3D case
as,

$\displaystyle P_{2s}^{2}$ $\displaystyle=P_{l}^{2}$ (13) $\displaystyle
P_{3s}^{2}$ $\displaystyle=\frac{3}{4}P_{l}^{2}+\frac{1}{4}P_{d}^{2}$ (14)

Thus, one can notice from equation (14) that for the 3D case, the proposed
degree of shearing constitute not only the degree of linear polarization,
$P_{l}$, but also the degree of directionality, $P_{d}$, which measure the
stability of the plane that contains the polarization ellipse, or
equivalently, the measure of closeness of the state represented by
$\mathbf{\Phi}$ to that of a 2D state [19]. In a similar way, this can be
further extended for $n>3$. However, one needs an appropriate physical
interpretation of such extension to higher dimensions.

Further, using the anti-symmetric matrix, $\mathbf{\Phi}_{2}$, we express the
degree of circular polarization exactly as defined by Gil et al., as,

$P_{c}=\frac{\sqrt{2}{\left\lVert\mathbf{\Phi}_{2}\right\rVert}_{\text{F}}}{\operatorname{tr}(\mathbf{\Phi})}.$
(15)

This quantity measures all contributions to circular polarization and is also
invariant under unitary transformation.

Finally, we express the generalized degree of polarimetric purity by combining
the degree of shearing (12) and the degree of circular polarization (15), as,

$P_{nD}=\sqrt{P_{ns}^{2}+\frac{n}{2\,(n-1)}P_{c}^{2}}.$ (16)

Furthermore, one can relate the generalized degree of polarimetric purity for
2D and 3D cases to the three CPP parameters using equations (13), and (14),
and equation (15) as,

$\displaystyle P_{2D}$ $\displaystyle=\sqrt{P_{2s}^{2}+P_{c}^{2}}$ (17)
$\displaystyle=\sqrt{P_{l}^{2}+P_{c}^{2}}$ (18)

and

$\displaystyle P_{3D}$ $\displaystyle=\sqrt{P_{3s}^{2}+\frac{3}{4}P_{c}^{2}}$
(19)
$\displaystyle=\sqrt{\frac{3}{4}P_{l}^{2}+\frac{1}{4}P_{d}^{2}+\frac{3}{4}P_{c}^{2}}.$
(20)

In [19, 20], Gil has explicitly shown the relationships of $P_{2D}$ and
$P_{3D}$ with the CPP parameters. Therefore, $P_{3s}$ provides fractional
contributions from both $P_{l}$ and $P_{d}$ (i.e., polarization fraction due
to mixed states), whereas $P_{2s}$ provides pure contribution from $P_{l}$.
Using the derivation proposed in Sheppard et al., [17], we can show that,

$\displaystyle P_{L}^{2}$
$\displaystyle=\frac{3}{4}P_{l}^{2}+\frac{1}{4}P_{d}^{2}-\frac{1}{4}P_{c}^{2},$
(21) $\displaystyle=P_{3s}^{2}-\frac{1}{4}P_{c}^{2}$ (22)

where $P_{L}$ is defined in [17] as the degree of total linear polarization,
i.e., the contribution from both the purely polarized and mixed state.

Now, expanding $P_{ns}$, and $P_{c}$ in the expression of $P_{nD}$ given in
equation (16) in terms of the Frobenius norm and the matrix trace, we find
that,

$\displaystyle P_{nD}^{2}$
$\displaystyle=\left(\dfrac{n}{n-1}\right)\Bigg{[}\frac{{\left\lVert\mathbf{\Phi}_{1}\right\rVert}^{2}_{\text{F}}+{\left\lVert\mathbf{\Phi}_{2}\right\rVert}^{2}_{\text{F}}}{\left(\operatorname{tr}{\mathbf{\Phi}}\right)^{2}}\Bigg{]}$
(23)
$\displaystyle=\left(\dfrac{n}{n-1}\right)\Bigg{[}\frac{ns^{2}}{n^{2}m^{2}}\Bigg{]},$
(24)

and therefore,

$P_{nD}=\frac{s}{\sqrt{n-1}\,m}.$ (25)

which is coincident with equation (5).

Hence, we suitably verify the equivalence among the two approaches for the
expression of the generalized degree of polarimetric purity, $P_{nD}$.

##  4 Acknowledgement

The authors would like to thank Dr. David Roberts, the University of Adelaide,
Australia, for the excellent blog article
(https://thehighergeometer.wordpress.com/2018/07/28/what-is-the-curl-of-a-
vector-field-really/, accessed on 09/06/2021) which prompted the initial
conception of Approach II of this work.

## References

  * Brosseau [1998] C. Brosseau, _Fundamentals of polarized light: a statistical optics approach_. Wiley-Interscience, 1998.
  * Samson [1973] J. C. Samson, “Descriptions of the Polarization States of Vector Processes: Applications to ULF Magnetic Fields,” _Geophysical Journal International_ , vol. 34, no. 4, pp. 403–419, 12 1973. [Online]. Available: https://doi.org/10.1111/j.1365-246X.1973.tb02404.x
  * Barakat [1977] R. Barakat, “Degree of polarization and the principal idempotents of the coherency matrix,” _Opt. Commun._ , vol. 23, no. 2, pp. 147–150, 1977.
  * Ellis and Dogariu [2005] J. Ellis and A. Dogariu, “On the degree of polarization of random electromagnetic fields,” _Optics communications_ , vol. 253, no. 4-6, pp. 257–265, 2005.
  * Réfrégier et al. [2006] P. Réfrégier, M. Roche, and F. Goudail, “Invariant polarimetric contrast parameters of light with gaussian fluctuations in three dimensions,” _JOSA A_ , vol. 23, no. 1, pp. 124–133, 2006.
  * Dennis [2007] M. R. Dennis, “A three-dimensional degree of polarization based on rayleigh scattering,” _JOSA A_ , vol. 24, no. 7, pp. 2065–2069, 2007.
  * Setälä et al. [2002] T. Setälä, A. Shevchenko, M. Kaivola, and A. T. Friberg, “Degree of polarization for optical near fields,” _Physical Review E_ , vol. 66, no. 1, p. 016615, 2002.
  * Luis [2005] A. Luis, “Degree of polarization for three-dimensional fields as a distance between correlation matrices,” _Optics communications_ , vol. 253, no. 1-3, pp. 10–14, 2005.
  * Gil et al. [2018] J. J. Gil, A. Norrman, A. T. Friberg, and T. Setälä, “Polarimetric purity and the concept of degree of polarization,” _Physical Review A_ , vol. 97, no. 2, p. 023838, 2018.
  * San José and Gil [2011] I. San José and J. J. Gil, “Invariant indices of polarimetric purity: generalized indices of purity for n$\times$ n covariance matrices,” _Optics communications_ , vol. 284, no. 1, pp. 38–47, 2011.
  * Gil et al. [2004] J. Gil, J. Correas, P. Melero, and C. Ferreira, “Generalized polarization algebra,” _Monog. Sem. Mat. G. Galdeano_ , vol. 31, pp. 161–167, 2004.
  * Gil [2007] J. J. Gil, “Polarimetric characterization of light and media,” _Eur. Phys. J. Appl. Phys_ , vol. 40, pp. 1–47, 2007.
  * Bhattacharya et al. [2022] A. Bhattacharya, S. Dey, and A. C. Frery, “Scattering purity and complexity in radar polarimetry,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2022.
  * Gilmore [2012] R. Gilmore, _Lie groups, Lie algebras, and some of their applications_. Courier Corporation, 2012\.
  * Wolkowicz and Styan [1980] H. Wolkowicz and G. P. H. Styan, “Bounds for eigenvalues using traces,” _Linear Algebra and Its Applications_ , vol. 29, pp. 471–506, 1980.
  * Romano and Price [2012] J. D. Romano and R. H. Price, “Why no shear in “div, grad, curl, and all that”?” _American Journal of Physics_ , vol. 80, no. 6, pp. 519–524, 2012\.
  * Sheppard et al. [2022] C. J. R. Sheppard, A. Bendandi, A. L. Gratiet, and A. Diaspro, “Purity of 3d polarization,” _J. Opt. Soc. Am. A_ , vol. 39, no. 1, pp. 6–16, Jan 2022\. [Online]. Available: http://www.osapublishing.org/josaa/abstract.cfm?URI=josaa-39-1-6
  * Gil [2014] J. J. Gil, “Interpretation of the coherency matrix for three-dimensional polarization states,” _Physical Review A_ , vol. 90, no. 4, p. 043858, 2014\.
  * Gil [2016] ——, “Components of purity of a three-dimensional polarization state,” _JOSA A_ , vol. 33, no. 1, pp. 40–43, 2016.
  * Gil and Ossikovski [2017] J. J. Gil and R. Ossikovski, _Polarized light and the Mueller matrix approach_. CRC press, 2017.