# Characterizing cylindrical particles upon local measurements of two Stokes
parameters

Jon Lasa-Alonso<EMAIL_ADDRESS>Donostia International Physics Center,
Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain.  Centro de
Física de Materiales, Paseo Manuel de Lardizabal 5, 20018 Donostia-San
Sebastian, Spain.  Iker Gómez-Viloria  Centro de Física de Materiales, Paseo
Manuel de Lardizabal 5, 20018 Donostia-San Sebastian, Spain.  Álvaro Nodar
Centro de Física de Materiales, Paseo Manuel de Lardizabal 5, 20018 Donostia-
San Sebastian, Spain.  Aitzol García-Etxarri  Donostia International Physics
Center, Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain.
IKERBASQUE, Basque Foundation for Science, María Díaz de Haro 3, 48013 Bilbao,
Spain.  Gabriel Molina-Terriza  Donostia International Physics Center, Paseo
Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain.  Centro de Física
de Materiales, Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastian,
Spain.  IKERBASQUE, Basque Foundation for Science, María Díaz de Haro 3, 48013
Bilbao, Spain.  Jorge Olmos-Trigo<EMAIL_ADDRESS>Centro de Física de
Materiales, Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastian, Spain.

###### Abstract

Researchers routinely characterize optical samples by computing the scattering
cross-section. However, the experimental determination of this magnitude
requires the measurement and integration of the components of the scattered
field in all directions. Here, we propose a method to determine the scattering
cross-section and global polarization state of radiation through measurements
of two Stokes parameters at an angle of choice in far-field. The method
applies to cylindrically symmetric samples whose optical response is well-
described by a single multipolar order $j$. Moreover, the formalism is
applicable for a wide range of different illuminations, and it only requires
the use of a single camera and conventional wave plates. Our findings
significantly reduce the complexity of routine characterization measurements
for cylindrical samples in optical laboratories.

The conservation of energy implies that it cannot be created or destroyed,
only transferred or transformed. This principle is essential to understand the
behavior of various physical systems, including the scattering of
electromagnetic radiation. In a general linear electromagnetic scattering
problem, energy is extinguished in two ways: scattering and absorption.
Extinction refers to the reduction in the intensity of a light beam as it
passes through a medium, scattering is the deflection of light in different
directions, and absorption is the process by which a sample consumes energy,
usually transforming it into heat. These processes are all governed by energy
conservation and can be described quantitatively using cross-sections. The
cross-sections can be understood as the area of a sample that is effectively
available for light to interact with. As a result, the cross-sections depend
on the sample but also on the incident illumination. The conservation of
energy is usually expressed in terms of the extinction, $\sigma_{\rm{ext}}$,
scattering, $\sigma_{\rm{sca}}$, and absorption cross-sections,
$\sigma_{\rm{abs}}$, of the sample through Bohren and Huffman (2008):

$\sigma_{\rm{ext}}=\sigma_{\rm{sca}}+\sigma_{\rm{abs}}.$ (1)

Thus, measuring the scattering and absorption cross-sections is crucial to
understand how energy transfers from light to matter. For example, the optical
theorem relates the extinction cross-section of a scatterer with the forward
scattering amplitude Novotny and Hecht (2012); Evlyukhin  _et al._ (2016).
Thus, the optical theorem greatly simplifies the measurement of
$\sigma_{\text{ext}}$ for samples under plane wave illumination. In
particular, for lossless samples, the optical theorem also gives access to the
scattering cross-section since in these situations
$\sigma_{\rm{ext}}=\sigma_{\rm{sca}}$. However, it is essential to note that
the optical theorem has some limitations: it neglects the complex structure of
focused beams, limiting its applicability to plane wave illumination Newton
(1976); Lock  _et al._ (1995); Krasavin  _et al._ (2018). On the other hand,
for lossy scatterers, direct measurements of the absorption or scattering
cross-sections are required to determine how light is extinguished in the
scattering process. Thus, in general, the characterization of linear
scatterers requires capturing the extinction, scattering, or absorption cross-
sections through alternative means.

One method to determine the absorption cross-section is measuring temperature
changes in the sample during the interaction with the incident beam (the Joule
effect) Baffou  _et al._ (2010); Giannini  _et al._ (2011). This method is
demanding as the absorption cross-section does not have a universal
relationship with temperature, and it is highly dependent on the material,
shape, and size of the scatterer Baffou and Quidant (2013). An alternative
approach is to capture the power scattered by a sample in order to determine
its scattering cross-section. However, measuring the scattering cross-section
is also challenging as it requires integrating all the scattered field
components. In any case, obtaining the extinction, scattering, and absorption
cross-sections of samples is essential for many applications spanning
different fields of photonics. We can find examples of the pivotal role of
scattering and absorption cross-sections in the context of resonant
nanoparticles Luk’yanchuk  _et al._ (2010); García-Etxarri  _et al._ (2011);
Kuznetsov  _et al._ (2012, 2016); Luk’yanchuk  _et al._ (2015), optical forces
Nieto-Vesperinas  _et al._ (2010); Gómez-Medina  _et al._ (2012a, b); Olmos-
Trigo  _et al._ (2019a) and torques Marston and Crichton (1984); Nieto-
Vesperinas (2015); Canaguier-Durand and Genet (2014); Xu and Nieto-Vesperinas
(2019), Kerker conditions Nieto-Vesperinas  _et al._ (2011); Geffrin  _et al._
(2012); Fu  _et al._ (2013); Person  _et al._ (2013); Liu  _et al._ (2014);
Olmos-Trigo  _et al._ (2020a); Lasa-Alonso  _et al._ (2023, 2022); Olmos-Trigo
_et al._ (2020b, c); Schmidt  _et al._ (2015); Xu  _et al._ (2020), optical
anapoles Miroshnichenko  _et al._ (2015); Luk’yanchuk  _et al._ (2017); Feng
_et al._ (2017); Ho and Scherrer (2013); Zanganeh  _et al._ (2021); Parker
_et al._ (2020); Wei  _et al._ (2016); Sanz-Fernández  _et al._ (2021), or
surface-enhanced circular dichroism spectroscopy Tang and Cohen (2010);
García-Etxarri and Dionne (2013); Ho  _et al._ (2017); Hu  _et al._ (2019);
Graf  _et al._ (2019); Solomon  _et al._ (2018); Lasa-Alonso  _et al._ (2020);
Olmos-Trigo  _et al._ (2023), among others.

In this work, we present a method to capture the scattering cross-section and
the global polarization state of light from local densities of the Stokes
parameters. That is, we show that by measuring both intensity $I(r,\theta)$
and degree of circular polarization $V(r,\theta)$ of the scattered field, at
any particular angle of choice $\theta$ in far-field, we can obtain the
scattering cross-section and global polarization state of light for
cylindrically symmetric samples. Importantly, our proposed method only
requires a single camera and a few wave plates placed at a fixed position far
from the scatterer. In principle, our approach applies when the electric and
magnetic optical response of the sample can be fully described by one
multipolar order $j$. However, we show that, under the appropriate
illumination conditions, all-optical resonances of generic cylindrical
particles can be determined through our proposed technique. Therefore, our
findings are relevant for any photonic branch dealing with the interaction
between electromagnetic waves and cylindrically symmetric samples Luk’yanchuk
_et al._ (2010); García-Etxarri  _et al._ (2011); Kuznetsov  _et al._ (2012,
2016); Luk’yanchuk  _et al._ (2015); Nieto-Vesperinas  _et al._ (2010); Gómez-
Medina  _et al._ (2012a); Olmos-Trigo  _et al._ (2019a); Gómez-Medina  _et
al._ (2012b); Marston and Crichton (1984); Nieto-Vesperinas (2015); Canaguier-
Durand and Genet (2014); Xu and Nieto-Vesperinas (2019); Nieto-Vesperinas  _et
al._ (2011); Geffrin  _et al._ (2012); Fu  _et al._ (2013); Person  _et al._
(2013); Liu  _et al._ (2014); Olmos-Trigo  _et al._ (2020a, b, c); Schmidt
_et al._ (2015); Xu  _et al._ (2020); Miroshnichenko  _et al._ (2015); Lasa-
Alonso  _et al._ (2022, 2023); Luk’yanchuk  _et al._ (2017); Feng  _et al._
(2017); Ho and Scherrer (2013); Zanganeh  _et al._ (2021); Parker  _et al._
(2020); Wei  _et al._ (2016); Sanz-Fernández  _et al._ (2021); Tang and Cohen
(2010); García-Etxarri and Dionne (2013); Ho  _et al._ (2017); Hu  _et al._
(2019); Graf  _et al._ (2019); Solomon  _et al._ (2018); Lasa-Alonso  _et al._
(2020); Olmos-Trigo  _et al._ (2023).

It is a well-known result from electrodynamics that the electromagnetic fields
radiated by an arbitrary linear scatterer can be expanded in terms of
multipolar fields. Following Jackson’s notation in its third edition Jackson
(1999), one can write:

$\displaystyle\mathbf{E}(\mathbf{r})=Z\sum_{j,m}\left[ia_{E}(j,m)\mathbf{N}^{h}_{jm}(\mathbf{r})+a_{M}(j,m)\mathbf{M}^{h}_{jm}(\mathbf{r})\right],$
$\displaystyle
iZ\mathbf{H}(\mathbf{r})=\frac{1}{k}\nabla\times\mathbf{E}(\mathbf{r}).$ (2)

Here $\mathbf{E}(\mathbf{r})$ and $\mathbf{H}(\mathbf{r})$ are the scattered
electric and magnetic fields, respectively; $Z=\sqrt{\mu/\epsilon}$ is the
impedance of the medium in which the scatterer is embedded, $\mu$ and
$\epsilon$ being the magnetic permeability and electric permittivity,
respectively, and $k$ is the modulus of the wavevector associated to the
radiation. Moreover, $a_{E}(j,m)$ and $a_{M}(j,m)$ are the electric and
magnetic coefficients;
$\mathbf{M}^{h}_{jm}(\mathbf{r})=h^{\scriptscriptstyle(1)}_{j}(kr)\mathbf{X}_{jm}(\theta,\varphi)$
and
$\mathbf{N}^{h}_{jm}(\mathbf{r})=k^{-1}\nabla\times\mathbf{M}^{h}_{jm}(\mathbf{r})$
are Hansen multipoles, with $h^{\scriptscriptstyle(1)}_{j}(kr)$ the spherical
Hankel functions of the first kind and $\mathbf{X}_{jm}(\theta,\varphi)$ the
vector spherical harmonics as defined in Ref. Jackson (1999). Here $j$ and $m$
denote the multipolar order and the eigenvalue of the total angular momentum
in the $z$ direction, respectively.

The description of the scattered fields given in Eq. (2) is general and
applies to an extensive range of different radiating systems. However, for the
aim of this work, we will focus on cylindrical samples illuminated by beams
with a well-defined total angular momentum $m$. In particular, we will tackle
scatterers that may be well-described by a single multipole order $j$ in Eq.
(2). This implies that the response of the sample is given by both electric
and magnetic scattering coefficients with $j=1$ (dipolar), $j=2$
(quadrupolar), $j=3$ (octupolar)… etc. On the other hand, when particularizing
to cylindrically symmetric samples illuminated by beams with well-defined
total angular momentum, we are also fixing the value of parameter $m$ in the
sum given by Eq. (2). The physical reason behind this phenomenon resides in
the fact that, because of their symmetry, cylindrical samples preserve the $z$
component of the total angular momentum upon light-matter interactions. Hence,
it is direct to check that when considering this specific type of scatterers,
we can get rid of the summation over both $j$ and $m$ indices in the
expression of the scattered field given by Eq. (2).

In addition, if we switch to the helicity basis, i.e.
$\mathbf{E}^{\pm}(\mathbf{r})=\left[\mathbf{E}(\mathbf{r})\pm
iZ\mathbf{H}(\mathbf{r})\right]/2$, we are led to the following compact
expression of the scattered fields for this particular type of scatterers:

$\mathbf{E}^{\pm}(\mathbf{r})=a_{\pm}(j,m)\left(\frac{\mathbf{N}^{h}_{jm}(\mathbf{r})\pm\mathbf{M}^{h}_{jm}(\mathbf{r})}{\sqrt{2}}\right).$
(3)

Here, $\mathbf{E}^{\pm}(\mathbf{r})$ are proportional to the monochromatic
Riemann- Silberstein vectors Bialynicki-Birula (1996); Bialynicki-Birula and
Bialynicka-Birula (2013) and $a_{\pm}(j,m)=Z[ia_{E}(j,m)\pm
a_{M}(j,m)]/\sqrt{2}$. From now on, the $(j,m)$ multipolar dependence will be
assumed in the scattering coefficients $a_{\pm}(j,m)$, namely, $a_{\pm}\equiv
a_{\pm}(j,m)$ . Our interest in choosing the helicity basis lies in the fact
that, in the far-field, the local fields are related to the Stokes parameters.
In this regard, fundamental magnitudes such as the total scattering cross-
section, proportional to the $I$ parameter, or the helicity expectation value,
constructed as the ratio $V/I$, can be obtained from the far-field expressions
of the scattered electromagnetic field in the helicity basis Crichton and
Marston (2000). In the far-field, we can write Eq. (3) as Carrascal  _et al._
(1991):

$\mathbf{E}^{\pm}_{\text{far}}(\mathbf{r})=-a_{\pm}\frac{e^{ikr}}{kr}\frac{(-i)^{j+1}}{\sqrt{j(j+1)}}\left(\frac{\bm{\xi}_{jm}(\theta,\varphi)\pm
i\bm{\eta}_{jm}(\theta,\varphi)}{\sqrt{2}}\right),$ (4)

where we have defined the complex vector functions
$\bm{\xi}_{jm}(\theta,\varphi)=r\nabla Y_{jm}(\theta,\varphi)$ and
$\bm{\eta}_{jm}(\theta,\varphi)=\hat{r}\times\bm{\xi}_{jm}(\theta,\varphi)$,
with $Y_{jm}(\theta,\varphi)$ the scalar spherical harmonics.

The measurable total scattered intensity $\langle I\rangle_{\rm{sca}}$ and
total degree of circular polarization $\langle V\rangle_{\rm{sca}}$ in the
far-field can be written as:

$\displaystyle\langle
I\rangle_{\rm{sca}}=(kr)^{2}\int~{}I_{\text{far}}(r,\theta)~{}d\Omega,$ (5)
$\displaystyle\langle
V\rangle_{\rm{sca}}=(kr)^{2}\int~{}V_{\text{far}}(r,\theta)~{}d\Omega,$ (6)

where
$I_{\text{far}}(r,\theta)=(|\mathbf{E}^{+}_{\text{far}}(\mathbf{r})|^{2}+|\mathbf{E}^{-}_{\text{far}}(\mathbf{r})|^{2})/|E_{0}|^{2}$
and
$V_{\text{far}}(r,\theta)=(|\mathbf{E}^{+}_{\text{far}}(\mathbf{r})|^{2}-|\mathbf{E}^{-}_{\text{far}}(\mathbf{r})|^{2})/|E_{0}|^{2}$
are, the local scattered intensity and the local degree of circular
polarization in the far-field, respectively. Note that $E_{0}$ is the
amplitude of the incident electric field and, as a result,
$I_{\text{far}}(r,\theta)$ and $V_{\text{far}}(r,\theta)$ are dimensionless.
For a sample whose response is well-described by a single multipolar order
$j$, the explicit computation of the integrals given in Eq. (5) and Eq. (6)
yields

$\displaystyle\langle
I\rangle_{\rm{sca}}=\frac{|a_{+}|^{2}+|a_{-}|^{2}}{|E_{0}|^{2}},$
$\displaystyle\langle
V\rangle_{\rm{sca}}=\frac{|a_{+}|^{2}-|a_{-}|^{2}}{|E_{0}|^{2}}.$ (7)

The results given in Eq. (7) have already been reported previously Jackson
(1999). However, their connection to the local intensity and local degree of
circular polarization densities for cylindrical samples has not been studied
in detail.

Let us explicitly compute the intensity and degree of circular polarization
densities in the far-field. For cylindrical samples well-described by a single
multipolar order $j$, $\langle I\rangle_{\rm{sca}}$ and $\langle
V\rangle_{\rm{sca}}$ are determined by the following fundamental relation (see
Appendix A):

$\displaystyle\begin{pmatrix}\langle I\rangle_{\rm{sca}}\\\ \langle
V\rangle_{\rm{sca}}\end{pmatrix}=\frac{(kr)^{2}j(j+1)}{f_{jm}^{2}(\theta)-g_{jm}^{2}(\theta)}\begin{pmatrix}f_{jm}(\theta)&g_{jm}(\theta)\\\
g_{jm}(\theta)&f_{jm}(\theta)\end{pmatrix}\begin{pmatrix}I_{\text{far}}(r,\theta)\\\
V_{\text{far}}(r,\theta)\end{pmatrix},$ (8)

with the real scalar functions
$f_{jm}(\theta)=|\bm{\xi}_{jm}(\theta,\varphi)|^{2}$ and
$g_{jm}(\theta)=\text{Im}[\bm{\xi}^{*}_{jm}(\theta,\varphi)\cdot\bm{\eta}_{jm}(\theta,\varphi)]$.
The expression given by Eq. (8) constitutes the main result of this work. It
indicates that, for cylindrical samples well-described by a single multipolar
order $j$, fundamental quantities such as the integrated scattered intensity
or total degree of circular polarization are related to their local densities
in the far-field computed at any particular scattering angle $\theta$. That
is, by measuring $I_{\text{far}}(r,\theta)$ and $V_{\text{far}}(r,\theta)$ in
one direction far from the illuminated sample, we can straightforwardly obtain
$\langle I\rangle_{\rm{sca}}$ and $\langle V\rangle_{\rm{sca}}$ by applying a
matrix that, in turn, can be easily evaluated once $j$ and $m$ are fixed
Jackson (1999). It is important to note that the matrix elements do not depend
on the optical response of the cylindrical scatterers (see Appendix B).

To get a deeper insight into the generality of Eq. (8), let us now introduce
one of the most illustrative examples in the nanophotonics community:
cylindrically symmetric dipolar objects under the illumination of a
circularly-polarized plane-wave Luk’yanchuk  _et al._ (2010); García-Etxarri
_et al._ (2011); Kuznetsov  _et al._ (2012, 2016); Luk’yanchuk  _et al._
(2015); Nieto-Vesperinas  _et al._ (2010); Gómez-Medina  _et al._ (2012a, b);
Olmos-Trigo  _et al._ (2019a); Marston and Crichton (1984); Nieto-Vesperinas
(2015); Canaguier-Durand and Genet (2014); Xu and Nieto-Vesperinas (2019);
Nieto-Vesperinas  _et al._ (2011); Geffrin  _et al._ (2012); Fu  _et al._
(2013); Person  _et al._ (2013); Liu  _et al._ (2014); Olmos-Trigo  _et al._
(2020a, c, b); Schmidt  _et al._ (2015); Xu  _et al._ (2020); Miroshnichenko
_et al._ (2015); Lasa-Alonso  _et al._ (2022, 2023); Luk’yanchuk  _et al._
(2017); Feng  _et al._ (2017); Ho and Scherrer (2013); Zanganeh  _et al._
(2021); Parker  _et al._ (2020); Wei  _et al._ (2016); Sanz-Fernández  _et
al._ (2021); Tang and Cohen (2010); García-Etxarri and Dionne (2013); Ho  _et
al._ (2017); Hu  _et al._ (2019); Graf  _et al._ (2019); Solomon  _et al._
(2018); Lasa-Alonso  _et al._ (2020); Olmos-Trigo  _et al._ (2023). This
particular physical system fixes $j=1$ (electric and magnetic dipolar
response) and $m=p$, where $p=\pm 1$ is the helicity of the incident beam. In
this scenario, it can be shown (see Appendix C) that Eq. (8) reduces to the
following expression:

$\displaystyle\begin{pmatrix}\langle I\rangle_{\text{sca}}\\\ \langle
V\rangle_{\text{sca}}\end{pmatrix}=\frac{16\pi(kr)^{2}}{3\sin^{4}\theta}\begin{pmatrix}1+\cos^{2}\theta&-2p\cos\theta\\\
-2p\cos\theta&1+\cos^{2}\theta\end{pmatrix}\begin{pmatrix}I_{\text{far}}(r,\theta)\\\
V_{\text{far}}(r,\theta)\end{pmatrix}.$ (9)

Equation (9) can be used as a roadmap to infer global properties of the
scattered field by cylindrical symmetry dipolar objects, such as $\langle
I\rangle_{\text{sca}}$ and $\langle V\rangle_{\text{sca}}$, from a single
measurement of local Stokes parameters. Then, from $\langle
I\rangle_{\text{sca}}$ and $\langle V\rangle_{\text{sca}}$, we can construct
widely employed magnitudes such as the scattering cross-section,
$\sigma_{\text{sca}}$, and the helicity expectation value,
$\langle\Lambda\rangle$ Lasa-Alonso  _et al._ (2022); Olmos-Trigo  _et al._
(2020b). More explicitly, these magnitudes can be fixed through the following
relations:

$\displaystyle k^{2}\sigma_{\text{sca}}=\langle I\rangle_{\text{sca}},$
$\displaystyle\langle\Lambda\rangle=\frac{\langle
V\rangle_{\text{sca}}}{\langle I\rangle_{\text{sca}}}.$ (10)

Equation (10) shows that one can compute the scattering cross-sections and the
helicity expectation value from local Stokes parameters, as we previously
anticipated. Note that the absorption cross-section, $\sigma_{\rm{abs}}$, can
be straightforwardly obtained if the extinction cross-section,
$\sigma_{\rm{ext}}$, is captured by other means (see Eq. (1)).

At this stage, we have all the necessary components to demonstrate the power
of our novel method. In this regard, we would like to highlight one of the
most significant contributions to the field of high-refractive-index (HRI),
“Magnetic Light” by Arseniy Kuznetsov et al. Kuznetsov  _et al._ (2012). In
this work, the authors showed the experimental dark-field scattering spectra
of spherical Silicon nanoparticles. Note that Silicon has losses in this
wavelength range. Furthermore, they presented the theoretical scattering and
extinction spectra calculated by Mie theory for spherical Silicon
nanoparticles of different sizes in free space within the visible spectral
range.

Figure 1: a) Scattering efficiency, $\sigma_{\text{sca}}/a^{2}$ (solid black),
and extinction efficiency, $\sigma_{\text{ext}}/a^{2}$ (dashed blue),
calculated from Mie theory for a Silicon nanosphere with radius $a=75$ nm in
the visible spectral range. The scattering cross-sections are calculated from
Eqs. (9)-(10) by computing $I_{\text{far}}(\theta,r)$ and
$V_{\text{far}}(\theta,r)$ at $\theta=80^{\circ}$ (dashed orange),
$\theta=90^{\circ}$ (dashed green) and $\theta=100^{\circ}$ (dashed red). b)
Absorption efficiency, $\sigma_{\text{abs}}/a^{2}$ (solid black), calculated
from Mie theory. Absorption efficiencies calculated from
$\sigma_{\text{ext}}-\sigma_{\text{sca}}(\theta)$ for $\theta=80^{\circ}$
(dashed orange), $\theta=90^{\circ}$ (dashed green), $\theta=100^{\circ}$
(dashed red). c) Helicity expectation value, $\langle\Lambda\rangle$,
calculated from Mie theory (solid black). The helicity expectation value at
different angles is calculated from Eqs. (9)-(10) by computing $I(\theta,r)$
and $V(\theta,r)$ at $\theta=80^{\circ}$ (dashed orange), $\theta=90^{\circ}$
(dashed green) and $\theta=100^{\circ}$ (dashed red).

In Fig. 1a-c, we depict several magnitudes for a Silicon nanosphere with
radius $a=75$ nm (Fig. 3c in Ref. Kuznetsov  _et al._ (2012)) under the
illumination of a circularly polarized plane-wave computed both from Mie
theory Mie (1908); Bohren and Huffman (2008) and through our novel method,
summarized in Eqs. (9)-(10) for dipolar objects. Let us just remark that the
results reported for different angles $\theta$ are associated with
measurements of the local Stokes parameters with a single camera placed far
from the sample. In Fig. 1a), we show the scattering efficiency,
$\sigma_{\text{sca}}/a^{2}$ (solid black), and the extinction efficiency,
$\sigma_{\text{ext}}/a^{2}$ (dashed blue), calculated from Mie theory.
Moreover, we show the scattering efficiency calculated from Eqs. (9)-(10) for
different scattering angles, namely, $\theta=80^{\circ}$ (dashed orange),
$\theta=90^{\circ}$ (dashed green) and $\theta=100^{\circ}$ (dashed red). The
agreement is total in the broadband wavelength interval of $475$ nm
$<\lambda<$ $800$ nm for every considered scattering angle. That is, when the
spherical particle can be fully described by electric and magnetic dipolar
modes. This agreement confirms that, indeed, global magnitudes such as the
scattering cross-section can be computed from local measurements of the $I$
and $V$ Stokes parameters.

In Fig. 1b), we show the absorption efficiency calculated from Mie theory and
from the relation
$\sigma_{\rm{abs}}(\theta)/a^{2}=\sigma_{\rm{ext}}/a^{2}-\sigma_{\rm{sca}}(\theta)/a^{2}$.
There is a reasonable agreement within the interval $475$ nm $<\lambda<$ $800$
nm. That is, when electric and magnetic modes can fully describe the
electromagnetic response of the spherical nanoparticle. This result indicates
that, under plane-wave illumination, one can obtain the absorption cross-
section without measuring temperature changes in spherical dipolar particles
during the interaction time. Note that this result holds as well for plasmonic
spherical dipolar particles.

Finally, in Fig. 1c), we depict the helicity expectation value computed from
both Mie theory and our novel approach. Note that we also capture
$\langle\Lambda\rangle$ with almost no error in the interval $475$ nm
$<\lambda<$ $800$ nm. Furthermore, we also show that it is possible to
identify the first Kerker condition Kerker  _et al._ (1983); Geffrin  _et al._
(2012); Fu  _et al._ (2013); Person  _et al._ (2013), at which helicity is
preserved Fernandez-Corbaton  _et al._ (2013); Olmos-Trigo  _et al._ (2019b);
Olmos-Trigo and Zambrana-Puyalto (2022) upon scattering, $\Lambda\sim 1$ at
$\lambda\sim 660$ nm.

Hitherto, we have shown that the scattering cross-section and the expected
value of helicity of dipolar particles under plane-wave illumination can be
determined through the relations specified in Eqs. (9)-(10). In this regard,
we have proved that the absorption cross-section can also be computed (via
extinction) from the scattering cross-section for these type of samples. In
other words, we have shown that our method permits the characterization of
small cylindrical particles under plane-wave illumination. In the forthcoming,
we show that our approach also allows the characterization of larger
cylindrically symmetric particles beyond the dipolar regime. For that aim, we
take advantage of the fact that the relation specified by Eq. (8) works under
a wider set of illumination conditions. In particular, we anticipate that
tightly-focused Laguerre-Gaussian (LG) beams with well-defined helicity, $p$,
and total angular momentum $m=\ell+p$, where $\ell$ is the orbital angular
momentum, can be employed for this task.

Let us now show all of this. In particular, we would like to drive attention
to another seminal contribution in the field of HRI nanoparticles: “Strong
magnetic response of submicron Silicon particles in the infrared” by García-
Etxarri et al García-Etxarri  _et al._ (2011). In this work, the authors
showed the scattering spectra calculated from Mie theory for a spherical
Silicon nanoparticle of radius $a=230$ nm in the near-infrared. Note that
Silicon has no losses in this wavelength range. We first excite such Silicon
nanosphere with a tightly-focused LG beam with $\ell=0$ and $p=1$ Gouesbet
_et al._ (1985, 2011). The focusing process is considered through an aplanatic
lens of numerical aperture $\rm{NA}=0.9$. As the total angular momentum for
this type of beam is $m=1$, we may still employ the relation given by Eq. (9).
Consequently, our method should facilitate the identification of dipolar
resonances. As we can infer from Fig. 2a), the scattering cross-section is
fully captured from local measurements at either $\theta=60^{\circ}$,
$\theta=90^{\circ}$ or $\theta=120^{\circ}$ in the wavelength interval of
$1300$ nm $<\lambda<2000$ nm. That is, when the scattering can be fully
described with electric and magnetic dipoles. However, due to the contribution
of the quadrupole, the scattering cross-section is ill-captured for
$\lambda<1300$ nm as the scattering cannot be fully described by a single
multipolar order $j$. Note that the disagreement in the local measurements at
$\lambda\sim 1300$ nm for different scattering angles is useful to identify
the regime at which the sample ceases to be dipolar.

Figure 2: Normalized scattering cross-section,
$\sigma_{\text{sca}}/\sigma_{\text{max}}$, for a Silicon nanosphere of radius
$a=230$ nm calculated from the Generalized Lorentz-Mie theory (solid black)
and at $\theta=60^{\circ}$ (dashed pink), $\theta=90^{\circ}$ (dashed green),
$\theta=120^{\circ}$ (dotted blue). $\lambda$ represents the wavelength of the
incident field. a) The incident field is a tightly-focused LG beam with total
angular momentum $m=1$ and helicity $p=1$. b) The incident field is a tightly-
focused LG beam with total angular momentum $m=2$ and helicity $p=1$. In both
cases, the numerical aperture is $\rm{NA}=0.9$.

At this point, as our method works for a wider range of illuminating fields,
we may change the incident beam and excite the scatterer with a tightly-
focused LG of $\ell=1$ and $p=1$ to fully capture the quadrupolar resonance.
The result is depicted in Fig. 2b). Such a remarkable outcome stems from the
fact that the illuminating beam has a total angular momentum $m=2$ and, thus,
it cannot excite any dipolar ($j=1$) resonances within the scatterer. It is a
well-known result that incident fields with total angular momentum $m$ cannot
excite multipolar resonances of order $j<|m|$ in cylindrically symmetric
scatterers Zambrana-Puyalto  _et al._ (2012). As a result, the underlying
dipolar resonances are cleared out, and our method remains valid for the
characterization of the particle in the range $\lambda<1300$ nm. Of course,
when exciting a quadrupolar mode ($j=2$) with an LG beam with total angular
momentum $m=2$, the matrix given by Eq. (8) has to be re-evaluated. In this
case, the functions that have to be employed are
$f_{22}(\theta)=(15/8\pi)\sin^{2}\theta(1+\cos^{2}\theta)$ and
$g_{22}(\theta)=-(15/4\pi)\sin^{2}\theta\cos\theta$. Note that this procedure
is recursive, i.e., an LG beam of $\ell=2$ and $p=1$ can be employed to
characterize octupolar ($j=3$) spectral regions and so on. Also, the numerical
aperture employed in Fig. 2a) and 2b) is the same. This implies that the
characterization method can be implemented within a single set up in which the
orbital angular momentum of the incident beam is switchable.

In conclusion, we have presented a novel method to capture the scattering
cross-section and global polarization state of radiation through local
measurements of the Stokes parameters at an angle of choice in the far-field.
This novel approach requires a single camera and conventional wave plates,
significantly reducing the complexity of routine characterization measurements
for cylindrically symmetric samples in optical laboratories. Moreover, the
method is applicable under a wide range of illumination conditions, which
implies that it can be used in a wide range of experimental settings. For
instance, we have shown its applicability for plane waves and LG beams of
different orbital angular momenta. In particular, we have shown that the use
of tightly-focused LG beams permits the spectral characterization of
multipolar cylindrical scatterers. Experimentally, the construction of LG
beams with different $\ell$ values can be carried out through the use of a
Spatial Light Modulator, which may easily be integrated in an optical set up.
Thus, our characterization method is completely feasible with the instruments
currently available.

Overall, our findings have important implications for the field of optics. By
simplifying and streamlining the characterization process, our method could
enable researchers to more accurately and efficiently study a wide range of
cylindrically symmetric samples in a plethora of photonic branches.

## Acknowledgments

J.O.T. acknowledges support from the Juan de la Cierva fellowship No.
FJC2021-047090-I of MCIN/AEI /10.13039/501100011033 and NextGenerationEU/PRTR.
J.O.T, J.L.A., and A.G.E acknowledge support from Project No.
PID2019-109905GA-C22 of the Spanish Ministerio de Ciencia, Innovación y
Universidades (MICIU). We received funding from the IKUR Strategy under the
collaboration agreement between Ikerbasque Foundation and DIPC on behalf of
the Department of Education of the Basque Government, as well the Basque
Government Elkartek program (KK-2021/00082) and the Centros Severo Ochoa
AEI/CEX2018-000867-S from the Spanish Ministerio de Ciencia e Innovación. A.N.
acknowledges the financial support from the Spanish Ministry of Science and
Innovation and the Spanish government agency of research
MCIN/AEI/10.13039/501100011033 through Project Ref. No. PID2019-107432GB-I00
and from the Department of Education, Research and Universities of the Basque
Government through Project Ref. No. IT1526-22. G.M.T received funding from the
IKUR Strategy under the collaboration agreement between Ikerbasque Foundation
and DIPC/MPC on behalf of the Department of Education of the Basque
Government.

##  Appendix A Response of a cylindrically symmetric object well-described by
a single multipolar order $j$

In this Appendix, we show how to derive the fundamental relation expressed in
Eq. (8) starting from the fields in Eq. (4). In the derivation we omit the
explicit dependence of the fields on the $(r,\theta,\phi)$ variables. By
taking the modulus squared of the fields expressed in Eq. (4), we arrive to:

$|\mathbf{E}_{\text{far}}^{\pm}|^{2}=|a_{\pm}|^{2}\left(\frac{|\bm{\xi}_{jm}|^{2}\mp\text{Im}\left[\bm{\xi}_{jm}^{*}\cdot\bm{\eta}_{jm}\right]}{(kr)^{2}j(j+1)}\right),$
(11)

where we have used that $|\bm{\eta}_{jm}|^{2}=|\bm{\xi}_{jm}|^{2}$.
Normalizing the expression by the modulus of the amplitude of the incident
electric field, $|E_{0}|^{2}$, and taking the sum and difference of the
helicity components expressed in Eq. (11) we arrive to:

$\displaystyle I_{\text{far}}$
$\displaystyle=\frac{1}{(kr^{2})j(j+1)}\left[f_{jm}\langle
I\rangle_{\text{sca}}-g_{jm}\langle V\rangle_{\text{sca}}\right]$ (12)
$\displaystyle V_{\text{far}}$
$\displaystyle=\frac{1}{(kr^{2})j(j+1)}\left[-g_{jm}\langle
I\rangle_{\text{sca}}+f_{jm}\langle V\rangle_{\text{sca}}\right],$ (13)

where
$I_{\text{far}}=(|\mathbf{E}^{+}_{\text{far}}|^{2}+|\mathbf{E}^{-}_{\text{far}}|^{2})/|E_{0}|^{2}$
and
$V_{\text{far}}=(|\mathbf{E}^{+}_{\text{far}}|^{2}-|\mathbf{E}^{-}_{\text{far}}|^{2})/|E_{0}|^{2}$
are the local Stokes parameters measured fram from the scatterer. Also, we
have defined the real scalar functions $f_{jm}=|\bm{\xi}_{jm}|^{2}$ and
$g_{jm}=\text{Im}[\bm{\xi}_{jm}^{*}\cdot\bm{\eta}_{jm}]$. Finally, we have
employed Eq. (7) of the main text to express the integrated Stokes parameters
$\langle I\rangle_{\text{sca}}$ and $\langle V\rangle_{\text{sca}}$.

Note that the relation expressed in Eqs. (12)-(13) can be compactly written in
matrix notation as:

$\begin{pmatrix}I_{\text{far}}\\\
V_{\text{far}}\end{pmatrix}=\frac{1}{(kr^{2})j(j+1)}\begin{pmatrix}f_{jm}&&-g_{jm}\\\
-g_{jm}&&f_{jm}\end{pmatrix}\begin{pmatrix}\langle I\rangle_{\text{sca}}\\\
\langle V\rangle_{\text{sca}}\end{pmatrix}.$ (14)

By taking the inverse of the relation expressed above we finally obtain:

$\begin{pmatrix}\langle I\rangle_{\text{sca}}\\\ \langle
V\rangle_{\text{sca}}\end{pmatrix}=\frac{(kr^{2})j(j+1)}{f_{jm}^{2}-g_{jm}^{2}}\begin{pmatrix}f_{jm}&&g_{jm}\\\
g_{jm}&&f_{jm}\end{pmatrix}\begin{pmatrix}I_{\text{far}}\\\
V_{\text{far}}\end{pmatrix},$ (15)

which is exactly the expression given in Eq. (8) of the main text.

##  Appendix B General form of $f_{jm}(\theta)$ and $g_{jm}(\theta)$

The general form of $f_{jm}(\theta)$ and $g_{jm}(\theta)$ functions is:

$f_{jm}(\theta)=C^{2}_{jm}\left\\{\sin^{2}\theta\left[\frac{\partial
P_{j}^{m}(\cos\theta)}{\partial(\cos\theta)}\right]^{2}+\frac{m^{2}}{\sin^{2}\theta}\left[P_{j}^{m}(\cos\theta)\right]^{2}\right\\}$
(16)

and

$g_{jm}(\theta)=2C_{jm}^{2}\left\\{m~{}P_{j}^{m}(\cos\theta)\frac{\partial
P_{j}^{m}(\cos\theta)}{\partial(\cos\theta)}\right\\},$ (17)

where $P_{j}^{m}(\cos\theta)$ are the associated Legendre polynomials. Also,
the coefficients $C_{jm}$ are defined as:

$C_{jm}=\sqrt{\frac{2j+1}{4\pi}\frac{(j-m)!}{(j+m)!}}.$ (18)

##  Appendix C Response of a cylindrically symmetric dipolar object

In this Appendix, we show the steps to arrive from Eq. (8) to Eq. (9). To
start with, let us explicitly write the functions $\bm{\xi}_{\ell
m}(\theta,\varphi)$ and $\bm{\eta}_{\ell m}(\theta,\varphi)$ for $\ell=1$. We
also fix that the total angular momentum is $m=p=\pm 1$. Following Jackson’s
book in its third edition Jackson (1999), we can write

$\displaystyle\bm{\xi}_{1m}(\theta,\varphi)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sqrt{\frac{3}{2\pi}}\left(\cos\theta
e^{im\varphi}\hat{{\bf u}}_{\theta}+ime^{im\varphi}\hat{{\bf
u}}_{\varphi}\right)$ (19) $\displaystyle\bm{\eta}_{1m}(\theta,\varphi)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sqrt{\frac{3}{2\pi}}\left(\cos\theta
e^{im\varphi}\hat{{\bf u}}_{\theta}-ime^{im\varphi}\hat{{\bf
u}}_{\varphi}\right)$ (20)

Now, let us compute $f_{1m}(\theta)=|\bm{\xi}_{1m}(\theta,\varphi)|^{2}$ and
$g_{1m}(\theta)=\text{Im}\left[\bm{\xi}^{*}_{1m}(\theta,\varphi)\cdot\bm{\eta}_{1m}(\theta,\varphi)\right]$.
After some algebra, we get

$\displaystyle f_{1m}(\theta)=\frac{3}{8\pi}\left(1+\cos^{2}\theta\right),$
$\displaystyle g_{1m}(\theta)=-\frac{3}{4\pi}m\cos\theta.$ (21)

Now, by inserting Eq. (21) into Eq. (8), it is straightforward to arrive to

$\displaystyle\begin{pmatrix}\langle I\rangle_{\text{sca}}\\\ \langle
V\rangle_{\text{sca}}\end{pmatrix}=\frac{16\pi(kr)^{2}}{3\sin^{4}\theta}\begin{pmatrix}1+\cos^{2}\theta&-2m\cos\theta\\\
-2m\cos\theta&1+\cos^{2}\theta\end{pmatrix}\begin{pmatrix}I_{\text{far}}(r,\theta)\\\
V_{\text{far}}(r,\theta)\end{pmatrix}$ (22)

Notice that in Eq. (22) we have made use of the trigonometrical identity
$(1+\cos^{2}\theta)^{2}-4\cos^{2}\theta=\sin^{4}\theta$. Finally, if we
substitute $m=p$, which is valid for both an incident circularly polarized
plane-wave and a circularly polarized focused Gaussian beam, we arrive at Eq.
(9).

## References

  * Bohren and Huffman (2008) C. F. Bohren and D. R. Huffman,  _Absorption and scattering of light by small particles_ (John Wiley & Sons, 2008).
  * Novotny and Hecht (2012) L. Novotny and B. Hecht,  _Principles of nano-optics_ (Cambridge university press, 2012).
  * Evlyukhin  _et al._ (2016) A. B. Evlyukhin, T. Fischer, C. Reinhardt, and B. N. Chichkov, Physical Review B 94, 205434 (2016).
  * Newton (1976) R. G. Newton, American Journal of Physics 44, 639 (1976).
  * Lock  _et al._ (1995) J. A. Lock, J. T. Hodges, and G. Gouesbet, JOSA A 12, 2708 (1995).
  * Krasavin  _et al._ (2018) A. V. Krasavin, P. Segovia, R. Dubrovka, N. Olivier, G. A. Wurtz, P. Ginzburg, and A. V. Zayats, Light: Science & Applications 7, 1 (2018).
  * Baffou  _et al._ (2010) G. Baffou, R. Quidant, and F. J. García de Abajo, ACS nano 4, 709 (2010).
  * Giannini  _et al._ (2011) V. Giannini, A. I. Fernández-Domínguez, S. C. Heck, and S. A. Maier, Chemical reviews 111, 3888 (2011).
  * Baffou and Quidant (2013) G. Baffou and R. Quidant, Laser & Photonics Reviews 7, 171 (2013).
  * Luk’yanchuk  _et al._ (2010) B. Luk’yanchuk, N. I. Zheludev, S. A. Maier, N. J. Halas, P. Nordlander, H. Giessen, and C. T. Chong, Nature materials 9, 707 (2010).
  * García-Etxarri  _et al._ (2011) A. García-Etxarri, R. Gómez-Medina, L. S. Froufe-Pérez, C. López, L. Chantada, F. Scheffold, J. Aizpurua, M. Nieto-Vesperinas, and J. J. Sáenz, Opt. Express 19, 4815 (2011).
  * Kuznetsov  _et al._ (2012) A. I. Kuznetsov, A. E. Miroshnichenko, Y. H. Fu, J. Zhang, and B. Luk’Yanchuk, Sci. Rep. 2, 492 (2012).
  * Kuznetsov  _et al._ (2016) A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk’yanchuk, Science 354, aag2472 (2016).
  * Luk’yanchuk  _et al._ (2015) B. S. Luk’yanchuk, N. V. Voshchinnikov, R. Paniagua-Domínguez, and A. I. Kuznetsov, ACS Photonics 2, 993 (2015).
  * Nieto-Vesperinas  _et al._ (2010) M. Nieto-Vesperinas, J. Sáenz, R. Gómez-Medina, and L. Chantada, Opt. Express 18, 11428 (2010).
  * Gómez-Medina  _et al._ (2012a) R. Gómez-Medina, B. García-Cámara, I. Suárez-Lacalle, L. Froufe-Pérez, F. González, F. Moreno, M. Nieto-Vesperinas, and J. Sáenz, Photonics and Nanostructures-Fundamentals and Applications 10, 345 (2012a).
  * Gómez-Medina  _et al._ (2012b) R. Gómez-Medina, L. Froufe-Pérez, M. Yépez, F. Scheffold, M. Nieto-Vesperinas, and J. J. Sáenz, Phys. Rev. A 85, 035802 (2012b).
  * Olmos-Trigo  _et al._ (2019a) J. Olmos-Trigo, C. Sanz-Fernández, F. S. Bergeret, and J. J. Sáenz, Opt. Lett. 44, 1762 (2019a).
  * Marston and Crichton (1984) P. L. Marston and J. H. Crichton, Phys. Rev. A 30, 2508 (1984).
  * Nieto-Vesperinas (2015) M. Nieto-Vesperinas, Phys. Rev. A 92, 043843 (2015).
  * Canaguier-Durand and Genet (2014) A. Canaguier-Durand and C. Genet, Physical Review A 89, 033841 (2014).
  * Xu and Nieto-Vesperinas (2019) X. Xu and M. Nieto-Vesperinas, Physical review letters 123, 233902 (2019).
  * Nieto-Vesperinas  _et al._ (2011) M. Nieto-Vesperinas, R. Gomez-Medina, and J. J. Saenz, J. Opt. Soc. Am. A 28, 54 (2011).
  * Geffrin  _et al._ (2012) J.-M. Geffrin, B. García-Cámara, R. Gómez-Medina, P. Albella, L. S. Froufe-Pérez, C. Eyraud, A. Litman, R. Vaillon, F. González, M. Nieto-Vesperinas, J. J. Sáenz, and F. Moreno, Nat. Commun. 3, 1171 (2012).
  * Fu  _et al._ (2013) Y. H. Fu, A. I. Kuznetsov, A. E. Miroshnichenko, Y. F. Yu, and B. Luk’yanchuk, Nat. Commun. 4, 1527 (2013).
  * Person  _et al._ (2013) S. Person, M. Jain, Z. Lapin, J. J. Sáenz, G. Wicks, and L. Novotny, Nano Lett. 13, 1806 (2013).
  * Liu  _et al._ (2014) W. Liu, J. Zhang, B. Lei, H. Ma, W. Xie, and H. Hu, Opt. Express 22, 16178 (2014).
  * Olmos-Trigo  _et al._ (2020a) J. Olmos-Trigo, C. Sanz-Fernández, D. R. Abujetas, J. Lasa-Alonso, N. de Sousa, A. García-Etxarri, J. A. Sánchez-Gil, G. Molina-Terriza, and J. J. Sáenz, Physical Review Letters 125, 073205 (2020a).
  * Lasa-Alonso  _et al._ (2023) J. Lasa-Alonso, J. Olmos-Trigo, C. Devescovi, P. Hernández, A. García-Etxarri, and G. Molina-Terriza, “Resonant helicity mixing of electromagnetic waves propagating through matter,” (2023), arXiv:2302.02647 [physics.optics] .
  * Lasa-Alonso  _et al._ (2022) J. Lasa-Alonso, J. Olmos-Trigo, A. García-Etxarri, and G. Molina-Terriza, Mater. Adv. 3, 4179 (2022).
  * Olmos-Trigo  _et al._ (2020b) J. Olmos-Trigo, D. R. Abujetas, C. Sanz-Fernández, X. Zambrana-Puyalto, N. de Sousa, J. A. Sánchez-Gil, and J. J. Sáenz, Physical Review Research 2, 043021 (2020b).
  * Olmos-Trigo  _et al._ (2020c) J. Olmos-Trigo, D. R. Abujetas, C. Sanz-Fernández, J. A. Sánchez-Gil, and J. J. Sáenz, Physical Review Research 2, 013225 (2020c).
  * Schmidt  _et al._ (2015) M. K. Schmidt, J. Aizpurua, X. Zambrana-Puyalto, X. Vidal, G. Molina-Terriza, and J. J. Sáenz, Phys. Rev. Lett. 114, 113902 (2015).
  * Xu  _et al._ (2020) X. Xu, M. Nieto-Vesperinas, C.-W. Qiu, X. Liu, D. Gao, Y. Zhang, and B. Li, Laser & Photonics Reviews 14, 1900265 (2020).
  * Miroshnichenko  _et al._ (2015) A. E. Miroshnichenko, A. B. Evlyukhin, Y. F. Yu, R. M. Bakker, A. Chipouline, A. I. Kuznetsov, B. Luk’yanchuk, B. N. Chichkov, and Y. S. Kivshar, Nature communications 6, 1 (2015).
  * Luk’yanchuk  _et al._ (2017) B. Luk’yanchuk, R. Paniagua-Domínguez, A. I. Kuznetsov, A. E. Miroshnichenko, and Y. S. Kivshar, Phys. Rev. A 95, 063820 (2017).
  * Feng  _et al._ (2017) T. Feng, Y. Xu, W. Zhang, and A. E. Miroshnichenko, Phys. Rev. Lett. 118, 173901 (2017).
  * Ho and Scherrer (2013) C. M. Ho and R. J. Scherrer, Physics Letters B 722, 341 (2013).
  * Zanganeh  _et al._ (2021) E. Zanganeh, A. Evlyukhin, A. Miroshnichenko, M. Song, E. Nenasheva, and P. Kapitanova, Physical Review Letters 127, 096804 (2021).
  * Parker  _et al._ (2020) J. A. Parker, H. Sugimoto, B. Coe, D. Eggena, M. Fujii, N. F. Scherer, S. K. Gray, and U. Manna, Physical Review Letters 124, 097402 (2020).
  * Wei  _et al._ (2016) L. Wei, Z. Xi, N. Bhattacharya, and H. P. Urbach, Optica 3, 799 (2016).
  * Sanz-Fernández  _et al._ (2021) C. Sanz-Fernández, M. Molezuelas-Ferreras, J. Lasa-Alonso, N. de Sousa, X. Zambrana-Puyalto, and J. Olmos-Trigo, Laser & Photonics Reviews , 2100035 (2021).
  * Tang and Cohen (2010) Y. Tang and A. E. Cohen, Physical review letters 104, 163901 (2010).
  * García-Etxarri and Dionne (2013) A. García-Etxarri and J. A. Dionne, Physical Review B 87, 235409 (2013).
  * Ho  _et al._ (2017) C.-S. Ho, A. Garcia-Etxarri, Y. Zhao, and J. Dionne, ACS Photonics 4, 197 (2017).
  * Hu  _et al._ (2019) J. Hu, M. Lawrence, and J. A. Dionne, ACS Photonics 7, 36 (2019).
  * Graf  _et al._ (2019) F. Graf, J. Feis, X. Garcia-Santiago, M. Wegener, C. Rockstuhl, and I. Fernandez-Corbaton, ACS Photonics 6, 482 (2019).
  * Solomon  _et al._ (2018) M. L. Solomon, J. Hu, M. Lawrence, A. García-Etxarri, and J. A. Dionne, ACS Photonics 6, 43 (2018).
  * Lasa-Alonso  _et al._ (2020) J. Lasa-Alonso, D. R. Abujetas, Á. Nodar, J. A. Dionne, J. J. Sáenz, G. Molina-Terriza, J. Aizpurua, and A. García-Etxarri, Acs Photonics 7, 2978 (2020).
  * Olmos-Trigo  _et al._ (2023) J. Olmos-Trigo, J. Lasa-Alonso, I. Gómez-Viloria, G. Molina-Terriza, and A. García-Etxarri, “Capturing near-field circular dichroism enhancements from far-field measurements,” (2023), arXiv:2301.01248 [physics.optics] .
  * Jackson (1999) J. D. Jackson,  _Classical Electrodynamics_ (John Wiley & Sons, New York, 1999).
  * Bialynicki-Birula (1996) I. Bialynicki-Birula,  _Progress in Optics_ , 36, 245 (1996).
  * Bialynicki-Birula and Bialynicka-Birula (2013) I. Bialynicki-Birula and Z. Bialynicka-Birula, Journal of Physics A: Mathematical and Theoretical 46, 053001 (2013).
  * Crichton and Marston (2000) J. H. Crichton and P. L. Marston, Electronic Journal of Differential Equations 4, 37 (2000).
  * Carrascal  _et al._ (1991) B. Carrascal, G. A. Estevez, P. Lee, and V. Lorenzo, European Journal of Physics 12, 184 (1991).
  * Mie (1908) G. Mie, Annalen der physik 330, 377 (1908).
  * Kerker  _et al._ (1983) M. Kerker, D.-S. Wang, and C. Giles, J. Opt. Soc. Am. A 73, 765 (1983).
  * Fernandez-Corbaton  _et al._ (2013) I. Fernandez-Corbaton, X. Zambrana-Puyalto, N. Tischler, X. Vidal, M. L. Juan, and G. Molina-Terriza, Phys. Rev. Lett. 111, 060401 (2013).
  * Olmos-Trigo  _et al._ (2019b) J. Olmos-Trigo, C. Sanz-Fernández, A. García-Etxarri, G. Molina-Terriza, F. S. Bergeret, and J. J. Sáenz, Phys. Rev. A 99, 013852 (2019b).
  * Olmos-Trigo and Zambrana-Puyalto (2022) J. Olmos-Trigo and X. Zambrana-Puyalto, Physical Review Applied 18, 044007 (2022).
  * Gouesbet  _et al._ (1985) G. Gouesbet, G. Gréhan, and B. Maheu, Journal of optics 16, 83 (1985).
  * Gouesbet  _et al._ (2011) G. Gouesbet, J. Lock, and G. Gréhan, Journal of Quantitative Spectroscopy and Radiative Transfer 112, 1 (2011).
  * Zambrana-Puyalto  _et al._ (2012) X. Zambrana-Puyalto, X. Vidal, and G. Molina-Terriza, Opt. Express 20, 24536 (2012).


# THE CALIFORNIA LEGACY SURVEY IV. LONELY, POOR, AND ECCENTRIC: A COMPARISON
BETWEEN SOLITARY AND NEIGHBORLY GAS GIANTS

Lee J. Rosenthal Cahill Center for Astronomy $\&$ Astrophysics, California
Institute of Technology, Pasadena, CA 91125, USA  Andrew W. Howard Cahill
Center for Astronomy $\&$ Astrophysics, California Institute of Technology,
Pasadena, CA 91125, USA  Heather A. Knutson Division of Geological and
Planetary Sciences, California Institute of Technology, Pasadena, CA 91125,
USA  Benjamin J. Fulton Cahill Center for Astronomy $\&$ Astrophysics,
California Institute of Technology, Pasadena, CA 91125, USA  IPAC-NASA
Exoplanet Science Institute, Pasadena, CA 91125, USA

(Received May 23rd, 2022; Revised Sep 12th, 2023; Accepted Oct 2nd, 2023)

###### Abstract

We compare systems with single giant planets to systems with multiple giant
planets using a catalog of planets from a high-precision radial velocity
survey of FGKM stars. Our comparison focuses on orbital properties, planet
masses, and host star properties. We use hierarchical methods to model the
orbital eccentricity distributions of giant singles and giant multis, and find
that the distributions are distinct. The multiple giant planets typically have
moderate eccentricities and their eccentricity distribution extends to
$e=0.47$ (90th percentile), while the single giant planets have a pile-up of
nearly circular orbits and a long tail that extends to $e=0.77$. We determine
that stellar hosts of multiple giants are distinctly more metal-rich than
hosts of solitary giants, with respective mean metallicities $0.228\pm 0.027$
vs. $0.129\pm 0.019$ dex. We measure the distinct occurrence distributions of
single and multiple giants with respect to orbital separation, and find that
single gas giants have a $\sim$2.3$\sigma$ significant hot ($a<0.06$) Jupiter
pile-up not seen among multi giant systems. We find that the median mass
($M\sin i$ ) of giants in multiples is nearly double that of single giants
(1.71 $M_{\textrm{J}}$ vs. 0.92 $M_{\textrm{J}}$ ). We find that giant planets
in the same system have correlated masses, analogous to the ‘peas in a pod’
effect seen among less massive planets.

exoplanets — catalogs

††journal: ApJS††software: All code used in this paper is available at
github.com/California-Planet-Search/rvsearch and
github.com/leerosenthalj/CLSIV. This research makes use of GNU Parallel
(Tange, 2011). We made use of the following publicly available Python modules:
astropy (Astropy Collaboration et al., 2013), matplotlib (Hunter, 2007),
numpy/scipy (van der Walt et al., 2011), pandas (McKinney, 2010), emcee
(Foreman-Mackey et al., 2013), and RadVel (Fulton et al., 2018).

##  1 Introduction

Giant planets are central players in planetary systems. Their dynamics affect
the presence of small inner planets (Zhu & Wu, 2018; Bryan et al., 2019;
Rosenthal et al., 2021a), and can reveal key information about the formation
and dynamical history of a system. We can also better understand the impact
that giant planets have on each other by studying their eccentricity
distributions, as planets with non-zero eccentricities must have reached their
states either by interacting with protoplanetary disks, or by undergoing
dynamical interactions with other planets or stellar companions (Dawson &
Murray-Clay, 2013; Petrovich & Tremaine, 2016). While secular interactions
between planets and planet-disk effects can account for giants excited to
moderate eccentricities ($0.2\leq e\leq 0.6$) (Dawson & Johnson, 2018), it is
more difficult to use them as an explanation for the highest-eccentricity
giants that have been discovered to date, which extend to $e>0.9$ (Blunt et
al., 2019).

Prior work has shown that highly eccentric giants can reach their states due
to strong gravitational scattering events, possibly with other giants
(Chatterjee et al., 2008). This would imply that giant multiplicity is a key
factor in understanding the dynamical evolution and final states of planetary
systems. Perhaps systems in which we currently only see one eccentric giant
began with multiple giant planets that scattered each other, one into a bound
eccentric orbit and one into engulfment or an unbound trajectory. We can
explore this theory by comparing the eccentricity distributions of systems
with one observed giant planet and multiple observed giant planets, as well as
their semi-major axis occurrence distributions and host star properties.
Recent theoretical work has provided giant planet hypotheses that are easily
testable with radial velocity (RV) surveys. For instance, Jackson et al.
(2021) used a synthetic population to show that if an RV search for outer
giant companions to warm Jupiters (10–200 days) finds that a substantial
fraction of these warm giants have outer companions, then we should seriously
consider a formation scenario in which secular interactions cause cold giants
to migrate inward and become stable warm Jupiters.

The California Legacy Survey (CLS, Rosenthal et al., 2021b) is well suited for
a comparative study of single and multiple giant planets. As a three-decade
long blind RV survey, it produced a sample that is appropriate for a variety
of occurrence measurements, and contains over a hundred giant ($M\sin i\
>0.1M_{J}$) planets, in both single and multiple configurations. In this
paper, we leverage this survey to compare and contrast lonely giants and multi
giant systems. In Section 2, we review the star and planet catalog of the CLS.
In Section 3, we describe our methods for computing planet occurrence,
hierarchically modeling eccentricity distributions, and comparing the host
star properties of distinct planetary samples. In Section 4, we present our
results. In Section 5, we discuss our findings and their context.

##  2 Survey Review

The California Legacy Survey (Rosenthal et al., 2021b) is a sample of RV-
observed FGKM stars and their associated planets, created in order to provide
a stellar and planetary catalog for occurrence studies. We approximated a
quantifiably complete survey by selecting HIRES-observed stars that were
originally chosen without bias towards a higher or lower than average
likelihood of hosting planets. That left us with 719 stars. We used an
iterative periodogram method to search for planet candidates, and performed
uniform vetting using magnetic activity indicators and instrument systematics
to identify false positives.

This left us with 178 planets in our sample, including 134 planets with $M\sin
i\ >0.1M_{J}$. 65 of these gas giants lack detected companions (‘lonely
giants’), and 69 belong to 31 multiple-giant systems. Figure 1 shows our giant
catalog in eccentricity and semi-major axis space, split between giant singles
and giant multis in one panel, and color-coded by host star metallicity in
another panel. Figure 2 shows the catalog in eccentricity and metallicity
space.

Figure 1: Left: Eccentricity with respect to semi-major axis for single giant
planets and multiple giant planets. Vertical bars show 15.9–84.1$\%$
confidence intervals. Circle radius is proportional to $M\sin i$. Right: Same
as left, but color-coded by host star metallicity. Curve shows estimate of
tidal disruption limit assuming that periastron at 0.03 au leads to disruption
(Dawson & Johnson, 2018). Figure 2: Planet eccentricity versus host star
metallicity for single and multiple giant systems.

We use $0.1M_{J}$ or 30 $M_{\oplus}$ as our cutoff for giant planets because
this is roughly the minimum mass at which runaway gas accretion can begin
(Lissauer et al., 2009). This is also near the threshold for planets to have
compositions dominated by gaseous envelopes (Lee, 2019; Thorngren et al.,
2016). These threshold bases set a more liberal definition of giant planet
than one motivated by dynamical influence on planet formation and migration,
but is more relevant to a planet’s growth history and captures sub-Saturn
planets that can still have an impact on smaller planets. Furthermore, our
sample is reasonably insensitive to our definition of this threshold. When we
lower the threshold from 30 $M_{\oplus}$ to 20 $M_{\oplus}$ the number of
single giant systems only decreases from 65 to 64, and the number of multi
giant systems only increases from 31 to 34.

##  3 Methods

###  3.1 Occurrence model

This work measures the occurrence of single and multiple giant planets as a
function of multiple orbital properties. Many studies of RV or transit surveys
use the intuitive occurrence measurement method known as ”inverse detection
efficiency” (Howard et al., 2012; Petigura et al., 2013). According to this
procedure, one estimates occurrence in a region of parameter space by counting
up the planets found in that region, with each planet weighted by the search
completeness in that region. One can measure the search completeness map of a
survey by injecting many synthetic signals into each dataset, and computing
the fraction of signals in a given region that are recovered by the search
algorithm in use. Inverse detection efficiency is actually a specific case of
a Poisson likelihood method, in which one models an observed planet catalog as
the product of an underlying Poisson process and empirical completeness map
(Foreman-Mackey et al., 2014).

In this work, as in previous papers in the CLS (Fulton et al., 2021; Rosenthal
et al., 2021a), we used the Poisson likelihood method to model the occurrence
of giant planets, taking measurement uncertainty into account. Given an
observed population of observed planets with orbital and $M\sin i$ posteriors
$\\{\bm{\omega}\\}$, and associated survey completeness map $Q(\bm{\omega})$,
and assuming that our observed planet catalog is generated by a set of
independent Poisson process draws, we can evaluate a Poisson likelihood for a
given occurrence model $\Gamma(\bm{\omega}|\bm{\theta})$, where
$\bm{{\theta}}$ is a vector of model parameters. The observed occurrence
$\hat{\Gamma}(\bm{\omega}|\bm{\theta})$ of planets in our survey can be
modeled as the product of the measured survey completeness and some underlying
occurrence model,

$\hat{\Gamma}(\bm{\omega}|\bm{\theta})=Q(\bm{\omega})\Gamma(\bm{\omega}|\bm{\theta}).\\\
$ (1)

The Poisson likelihood for an observed population of objects is

$\mathcal{L}=e^{-\int\hat{\Gamma}(\bm{\omega}|\bm{\theta})\,d\bm{\omega}}\prod_{k=1}^{K}\hat{\Gamma}(\bm{\omega}_{k}|\bm{\theta}),\\\
$ (2)

where $K$ is the number of observed objects, and $\bm{\omega}_{k}$ is the
$k$th planet’s orbital parameter vector. The Poisson likelihood can be
understood as the product of the probability of detecting an observed set of
objects (the product term in Equation 2), and the probability of observing no
additional objects in the considered parameter space (the integral over
parameter space). Another way to understand this model is that it represents a
Poisson process, with some probability of generating planets with a given mass
and semi-major axis rate density.

Equations 1 and 2 serve as the foundation for our occurrence model, but do not
take into account uncertainty in our measurements of planetary orbits and
minimum masses. In order to do this, we use Markov Chain Monte Carlo methods
to empirically sample the orbital posteriors of each system (Rosenthal et al.,
2021b). We can hierarchically model the orbital posteriors of each planet in
our catalog by summing our occurrence model over many posterior samples for
each planet. The hierarchical Poisson likelihood is therefore approximated as

$\mathcal{L}\approx
e^{-\int\hat{\Gamma}(\bm{\omega}|\bm{\theta})\,d\bm{\omega}}\prod_{k=1}^{K}\frac{1}{N_{k}}\sum_{n=1}^{N_{k}}\frac{\hat{\Gamma}(\bm{\omega}_{k}^{n}|\bm{\theta})}{p(\bm{\omega}_{k}^{n}|\bm{\alpha})},\\\
$ (3)

where $N_{k}$ is the number of posterior samples for the $k$th planet in our
survey, and $\bm{\omega}_{k}^{n}$ is the $n$th sample of the $k$th planet’s
posterior. $p(\bm{\omega}|\bm{\alpha})$ is our prior on the individual planet
posteriors. We placed uniform priors on ln($M$sin$i$) and ln($a$). We used
emcee (Foreman-Mackey et al., 2013) to sample our hierarchical Poisson
likelihood.

###  3.2 Application to eccentricity

As done in prior eccentricity studies (Kipping, 2013; Van Eylen et al., 2019;
Bowler et al., 2020), we model the eccentricity distribution of a population
of exoplanets with the Beta distribution, because it is [0, 1] bound, flexible
in its shape, and has only two model parameters. The Beta distribution follows
the probability density function

$p(x|\alpha,\beta)\ =\
\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1}\\\
$ (4)

where $\Gamma(x)$ is the Gamma function, rather than the occurrence rate in
Equation 3. We only consider planets with $M\sin i$ $\geq$ 0.1
$M_{\textrm{J}}$ . In other words, we are only considering the eccentricity
distributions of planets with masses greater than a third of Saturn’s mass.

When we apply our hierarchical modeling framework to the eccentricity
posteriors of our giant planets, we produce the likelihood

$\mathcal{L}(\bm{e}|\bm{\theta})\propto\prod_{k=1}^{K}\frac{1}{N}\sum_{n=1}^{N}p(e_{k}^{n}|\bm{\theta}),\\\
$ (5)

where $\theta$ is the set of model parameters, $K$ is the number of planets,
$N$ is the number of samples drawn from each planet’s posterior, and
$e_{k}^{n}$ is the eccentricity of the $n$th sample drawn from the $k$th
planet’s posterior. The normalization term in Equation (3) disappears because
Equation (4) is a normalized probability density function.

###  3.3 Empirical distribution comparisons

In order to compare the host-star metallicities and other properties of
distinct groups of planetary systems, we use simple statistical tests to
compare the cumulative distribution functions, otherwise known as empirical
distributions, of these groups. Specifically, we use the Anderson-Darling
test, a statistical method for computing the probability that two empirical
distributions were drawn from the same underlying distribution.

##  4 Results

###  4.1 Single and multiple giants have distinct eccentricity distributions

We use the hierarchical model described in Section 3 to characterize the
distinct eccentricity distributions of single and multiple giant planets.
Figures 3 and 4 show our results when we fit a beta distribution to these
populations. For single giants, our model returns beta parameters
$\alpha=0.60^{+0.14}_{-0.12}$ and $\beta=1.45^{+0.34}_{-0.30}$. For multiple
giants, our model returns $\alpha=1.06^{+0.21}_{-0.18}$ and
$\beta=3.68^{+0.81}_{-0.70}$. While we define multiplicity as the number of
giant planets at any semi-major axis, our joint fits only include planets
beyond 0.3 au, since giants within that distance have short tidal
circularization timescales (Ogilvie, 2014). For reference, our population fit
to all giant planets beyond 0.3 au returns $\alpha=0.85^{+0.13}_{-0.11}$ and
$\beta=2.27^{+0.36}_{-0.32}$.

Both the beta distribution models and the raw posteriors show that single
giant planets have a pile-up of circular orbits and a long tail that extends
to $e\leq 1$, while the multiple giant planets have a more extended range of
moderate eccentricity and a sharp cutoff around $e\sim 0.7$. Hints of this
behavior are visible in in Figure 13 of Wright et al. (2009), which compared
the eccentricity distributions of all single and multi-planet systems
discovered to date and found a long single-planet tail not shared by the
multis. Bryan et al. (2016) reported a similar result with a larger planet
sample.

One possible explanation is that multiple giant systems have some probability
of experiencing secular interactions or scattering events, whereas truly
single giants need strong disk interactions or Kozai-Lidov interactions with a
stellar binary to reach high eccentricities. This would imply that the
observed population of more circular single giants have not experienced
companion-driven excitation, while the moderately eccentric multis and highly
eccentric singles represent different outcomes from planet-planet
interactions. Eccentric giants that appear to be lonely now may have been born
with nearby companions and scattered them into ejection or tidal engulfment,
while systems with two observed, moderately eccentric giants experienced
dynamical interactions that were strong enough to excite but not eject or
cause engulfment. The alternative of Kozai-Lidov excitation by binary stars is
entirely possible, particularly since a Gaia query reveals that most of the 12
CLS singles with e $>$ 0.5 have wide stellar companions.

Figure 3: Eccentricity distributions of single and multiple giants beyond 0.3
au. Purple/top distribution is for multiple giants, green/bottom distribution
is for single giants. Histograms are distributions of individual maximum
posterior values. Vertical bars show 68$\%$ confidence intervals

. Curves show the median and many draws from beta distribution posteriors,
produced using hierarchical inference.

Figure 4: Beta distribution posteriors for the single and multiple giant
eccentricity distributions. The two populations are $\sim 3\sigma$-distinct.

###  4.2 Hosts of multiple giants are more metal-rich than single giant hosts

We compare the host star metallicity distributions of systems with only one
observed giant planet and multiple observed giant planets to determine whether
these two groups are different. Figure 5 shows that hosts of multiple giant
planets are distinctly more metal rich than hosts of only one detected giant
planet, with 82.2$\%$ significance, that being the fraction of Anderson-
Darling tests performed on many samples from our distributions that fall
within a p-value $p\leq 0.01$. This finding is in agreement with the trend
first observed in Fischer & Valenti (2005). We split our giant host sample
along its median metallicity and found that the search completeness contours
for the two subsamples are nearly identical, as seen in Figure 6. This
provides confidence that the difference in these two distributions is not
caused by more extensive RV observations of our most metal-rich host stars.

We can interpret this phenomenon in several ways. One possible explanation is
that the probability of more than one giant planet forming around a star
increases with increasing host metallicity. This is consistent with
expectations from the core accretion theory of giant formation (Lissauer et
al., 2009). Another possible explanation is that metal-richer stars are more
likely to form multiple giant planets closer in, and that our sample of single
giants contains additional undetected giant companions beyond 30 au. This is
disfavored by the results from Fulton et al. (2021), which found evidence that
giant planets are less common beyond 8 au than within 1–8 au.

In light of our single-multi metallicity result, we explore whether the
eccentricity distribution or semi-major axis distribution of the giant planet
population is also metallicity-dependent. The hosts of circular ($e<0.2$)
giant planets and hosts of eccentric ($e\geq 0.2$) giant planets have median
metallicities separated by less than 1$\sigma$ in their standard errors, and
Anderson-Darling tests produce p-values almost entirely above 0.01. The same
is true of hot ($a<0.1$) giant hosts with respect to colder ($a\geq 0.1$)
giant hosts. This rules out the possibility that our single-multi metallicity
result is driven by $a$ or $e$ as confounding factors, rather than
multiplicity itself.

Splitting single giants along $e\geq 0.5$, further along the long tail of
eccentric singles, leads to inconclusive results. Figure 7 implies that there
is not a conclusive difference in host metallicity between single $e\geq 0.5$
hosts and hosts of less eccentric singles, and Anderson-Darling tests also
show no conclusive difference. The medians of the two groups are separated by
just 1.1$\sigma$. However, the same is all true of the difference between
single $e\geq 0.5$ hosts and all multi giant hosts; Anderson-Darling tests are
inconclusive, and the medians of the two groups are separated by 1$\sigma$. We
therefore conclude that our sample of single eccentric planets is not large
enough to differentiate its metallicity distribution from that of either the
population of single planets on circular orbits or the multi giant population.
This is not surprising, as our sample only contains 12 single planets with
eccentricities higher than 0.5. The two comparison samples are proportionally
larger, with 52 systems with single giants on low eccentricity orbits ($e$
$\leq$ 0.5) and 31 systems with multiple giant planets.

Figure 5: Cumulative metallicity distributions for hosts of single giant
planets and hosts of multiple giant planets. Solid and dashed steps show
distributions of median metallicity measurements, while transparent steps show
metallicities drawn many times from Gaussian distributions, with means and
standard deviations taken from measurement means and uncertainties. Anderson-
Darling tests show that the two groups are distinct.

Figure 6: Top: Average search completeness contours for giant hosts below
their median metallicity. Bottom: Average search completeness contours for
giant hosts above their median metallicity. The two contours are nearly
equivalent, implying that sensitivity bias to metal-rich stars is not the
cause of the single-multi host distinction seen in Figure 5. Figure 7:
Cumulative metallicity distributions for hosts of circular ($e<0.5$) giant
planets, hosts of eccentric ($e\geq 0.5$) giant planets, and multi giant
hosts. See Figure 5 for description. The eccentric single hosts are
inconclusively drawn from the same underlying metallicity distribution as the
less eccentric singles or the multis.

###  4.3 Hot Jupiter pile-up is missing from multis

We use our Poisson likelihood model to compute occurrence as a function of
orbital distance for the distinct single giant and multiple giant populations.
Figure 8 shows our results. There is a hot Jupiter pile-up within 0.6 au for
lonely giants, but no such pile-up for neighborly giants. We further
investigate this phenomenon by generating occurrence posteriors for three
classes of giants: warm Jupiters (0.1–0.4 au), moderately hot Jupiters
(0.06–0.1 au), and very hot Jupiters (0.03–0.06 au). The single giant and
multiple giant distributions of warm and hot Jupiters are indistinguishable,
but the very hot Jupiter distributions are separated with 97.6$\%$ confidence.
We measured this confidence as the fraction of random draws from the
difference between distributions that is greater than zero. When we define
multiplicity to include all massive companions, including brown dwarfs and
stellar binaries, one very hot Jupiter changes in multiplicity status, and our
separation confidence decreases to 94.4$\%$ significance.

Figure 8: Occurrence with respect to semi-major axis for single giant planets
and multiple giant planets. Steps and dots are median values, while vertical
lines are 15.9–84.1$\%$ confidence intervals. The single giant planets display
a pile-up of very hot Jupiters within 0.06 au, tentatively not shared by the
multiple giant planets.

While very hot Jupiters may be lonely, warm Jupiters ($\geq$ 0.1 au) may tell
another story. Jackson et al. (2021) predicts that warm Jupiters have outer
giant companions with significantly non-zero frequency, which could imply that
these companions interact secularly to bring the warm Jupiters to their final
states. Alternatively, Huang et al. (2016) found that warm Jupiters are more
likely to have small planet companions than hot Jupiters. This could imply
that warmer giants are less likely to have migrated inward, since migration
would disrupt the orbits of small companions. We test these hypotheses with
our sample of warm Jupiters, shown in Figure 9. In a sample of 21 systems, 10
show detected outer giant companions, while 11 do not. We satisfy the
prediction of Jackson et al. (2021), and therefore cannot rule out secular
migration. Further resolving this question requires a sample that is uniformly
sensitive to both cold gas giants and small rocky planets, which is not the
case for the CLS.

Figure 9: Lonely and neighborly warm Jupiters, described in Subsection 4.3. 11
systems have lonely warm Jupiters, while 10 systems host both a warm Jupiter
and at least one detected outer giant companion. This does not refute the
synthetically driven claim by Jackson et al. (2021) of secular warm Jupiter
migration.

###  4.4 Inconclusive mass functions

By investigating occurrence as a function of mass for single and multiple
giant planets separately, we may be able to infer whether the two groups form
via different pathways. Figure 10 shows occurrence rates with respect to
$M\sin i$ of both single and multiple giants within 0.1–2 au (multis can
include planets outside this range). These two distributions are $\leq
2\sigma$ separated in all four considered ranges of $M\sin i$ , but the multi
distribution tends toward higher masses than the single distribution. Leaving
aside our occurrence model, the median $M\sin i$ of our single giant planets
is 0.92 $M_{\textrm{J}}$, while the median $M\sin i$ of giant multiples is
1.71 $M_{\textrm{J}}$. While we cannot assume that the giant masses are
normally distributed, we can report that the standard error in the mean $M\sin
i$ for singles is 0.33 $M_{\textrm{J}}$, and that the standard error in the
mean $M\sin i$ for multis is 0.42 $M_{\textrm{J}}$.

Stepping back from our occurrence model, we can compare the cumulative
distributions of the two groups and make assumptions that the average
completeness contours for the single hosts and multi hosts are nearly
identical, similar to the metallicity-split contours in Figure 6. Figure 11
shows these cumulative distributions. AD tests show that there is a $>99.7\%$
chance of a p-value below 0.01 for matching single and multi giants within 3
au. This can be interpreted as the degree of confidence in the claim that
giants in multiples are consistently more massive than single giants.
Additionally, while eccentric singles appear to be more massive than circular
singles, this could be a selection effect of more eccentric detections with
greater mass.

Figure 10: Occurrence as a function of $M\sin i$ for single and multiple
giants. Steps and dots are median values, while vertical lines are
15.9–84.1$\%$ confidence intervals. The two distributions are not $\leq
2\sigma$ separated in all four considered ranges. Figure 11: Cumulative $M\sin
i$ distributions for multiple giants within 3 au and single giants split
across an eccentricity of 0.2. While the single and circular giants are
$>4\sigma$ separated from the multis, eccentric giants are inconclusively more
massive. Note that this may be a selection effect for more eccentric
detections with greater mass.

###  4.5 The properties of multiple giant planets

Comparative work requires taking a broader view of giant multiplicity, beyond
hot or warm Jupiters. Figure 12 shows the observed giant multiplicity
distribution of the CLS. Approximately two thirds of all stars that host a
giant planet host only one observed giant planet, while one third of giant
hosts show two or more observed giants. This distribution is not corrected for
search completeness; it only shows the multiplicity of detected giant planets.
This means that at least one third of stars that host one giant planet also
host two or more giant planets, unless we are somehow less sensitive to single
giants than to multi giants.

Figure 12: Observed multiplicity of giant planets in the CLS sample. Half as
many systems contain two or more giants as contain only one giant.

By looking at systems of multiple giant planets, we can also explore how
important global disk properties are for setting final planet masses compared
to the timeline and location of each individual planet’s formation. Focusing
on multiple giant systems, we investigate whether paired giant planets have
correlated masses. Specifically, we compare $M\sin i$ values of neighboring
giant planets, for each neighboring giant pair in a system. We measure a
Pearson correlation coefficient of $R=0.476\pm 0.010$ (error come from
bootstrapping measured $M\sin i$ values with uncertainties). We use two
independent methods to investigate whether sampling a population of
uncorrelated giant pairs may produce this value of $R$.

In one case, we use bootstrapping to approximate uncorrelated populations
using the observed CLS catalog. First, we repeatedly draw two giant planets
from the CLS catalog until we have thirty-one pairs, the same size as the
number of observed multi giant systems. We sample without replacement to
ensure that we are shuffling our entire giant catalog, rather than drawing the
same planets multiple times. Then, we measure the Pearson correlation
coefficient of this shuffled population. We also measure the standard
deviation in the log-ratio of planet pair $M\sin i$s, and assume coplanarity
to simplify this to $\sigma_{\mathrm{ln}(M_{2}/M_{1})}$. It must be noted that
coplanarity between giant planets is not a given, and that this limits the
significance of our results.

We repeat this process 105 times, until we build up bootstrapped distributions
of $R$ and $\sigma_{\mathrm{ln}(M_{2}/M_{1})}$. Our measured $R$ falls just
outside the 99.7$\%$ confidence interval of the bootstrapped distribution, at
the 99.71$\%$ confidence interval. This means that there is a 0.29$\%$
probability that an uncorrelated giant planet population could produce a
Pearson correlation coefficient of 0.45. Our measured
$\sigma_{\mathrm{ln}(M_{2}/M_{1})}$ falls on the 1.4$\%$ confidence interval
of the bootstrapped distribution.

We also use a synthetic population experiment to independently measure the
significance of our observed correlation. We drew giant planet pairs from a
uniform ln($M\sin i$ ) distribution and the semi-major axis distribution
reported in Fulton et al. (2021), with a minimum period semi-major axis ratio
ln($a_{2}/a_{1}$) $>$ 0.5. We impose this limit on the semi-major axis ratio
in order to prevent unlikely or unphysical draws of pairs within 3:2
resonance. We took our search completeness contours into account to only
select simulated planet pairs for which both planets are detected. We detail
our synthetic population algorithm in the list below.

  * •

Draw two values of ln($M\sin i$ ) from a uniform distribution between ln(30
$M_{\oplus}$)–ln(6000 $M_{\oplus}$), and two values of ln($a$) from the
nonparametric giant planet model reported in Fulton et al. (2021), bounded
between 0.1–10 au. We approximate this distribution as piecewise-uniform in
ln($a$), with a break at 1 au, integrated probability of 2/11 between 0.1–1
au, and integrated probability of 9/11 between 1–10 au.

  * •

Measure the search completeness and draw a random number from a uniform 0–1
distribution for each set of ln($a$) and ln($M\sin i$). For each set, if the
random number drawn is less than the detection probability at those
parameters, report the planet corresponding to that set as detected.

  * •

Keep the planet pair if we detect both planets and ln($a_{2}/a_{1}$) $>$ 0.5,
otherwise reject it.

  * •

Repeat until we have drawn the desired number of simulated planet pairs for
one synthetic population.

  * •

Repeat until we have drawn the desired number of synthetic populations.

With $10^{4}$ synthetic populations of 31 observed planet pairs each, this
test finds that our observed correlation falls just within the 99$\%$
confidence interval of the distribution that our uncorrelated simulations
produces.

##  5 Discussion

###  5.1 Metallicity, eccentricity, and multiplicity

Dawson & Murray-Clay (2013) found that giants with metal-rich hosts show
signatures of planet-planet scattering. In particular, giant planets between
0.1 AU and 1 AU are likely to be more eccentric when they orbit stars with
[Fe/H] $>$ 0\. This higher eccentricity could be caused by planet-planet
scattering events, which would imply that metal-rich stars are more likely to
host multiple giant planets. Our results support this hypothesis, as we find
that systems of multiple giant planets are more common around metal-rich
stars, and that systems with multiple giant planets have higher average
eccentricities than those in single giant planet systems. Table
LABEL:tab:comparisons shows the above properties side-by-side for the two
populations.

Table 1: Summary of population comparison between single giant systems and multi giant systems. Fraction of bound binaries comes from a literature estimate with Simbad and the Washington Double-star Catalog. Uncertainties come from Poisson counting error. Property | Characteristic | Single Giant Planets | Multiple Giant Planets
---|---|---|---
Eccentricity |  |  |
| 10th percentile, $e_{10\%}$ | 0.03 | 0.03
| Median, $e_{50\%}$ | 0.20 | 0.23
| 90th percentile, $e_{90\%}$ | 0.77 | 0.47
| Beta distribution $\alpha$ | $0.60^{+0.14}_{-0.12}$ | $1.06^{+0.21}_{-0.18}$
| Beta distribution $\beta$ | $1.45^{+0.34}_{-0.30}$ | $3.68^{+0.81}_{-0.70}$
Semi-major Axis |  |  |
| Comparison | HJ pile-up | HJ deficit
Mass |  |  |
| Median | 0.92 $M_{\textrm{J}}$ | 1.71 $M_{\textrm{J}}$
| Mean | 2.07 $\pm$ 0.33 $M_{\textrm{J}}$ | 2.88 $\pm$ 0.42 $M_{\textrm{J}}$
| Feature | Less massive | Intra-system uniformity
Iron Abundance |  |  |
|  $<[$Fe/H$]>$ |  $0.129\pm 0.019$ dex |  $0.228\pm 0.027$ dex
| Comparison | Enriched over Solar | More enriched
Stellar Binarity |  |  |
| Fraction Bound (Lit.) | $26.2\pm 6.3\%$ | $25.8\pm 9.1\%$

Although this story is clear, our data are less definitive on the question of
how the sub-population of single giants on highly eccentric orbits acquired
their eccentricities. There are few other mechanisms that could explain the
metal-poor, high-eccentricity tail of lonely giants above $e\sim 0.6$. The
most likely alternative to scattering is the Kozai-Lidov mechanism, by which a
massive outer companion secularly perturbs the orbit of an inner companion,
can potentially excite giants to $e>0.6$, but only given a high mutual
inclination and absence of other dynamical factors Naoz (2016). Even if this
could account for all of our single, highly eccentric giants, it would require
the presence of stellar binary companions or undetected giant companions in
these systems. We went back and examined the radial velocity data sets for the
eccentric single planets in our sample, and found that 4 out of 65 single
giant systems have $\geq 3\sigma$-significant parabolic or linear trends, and
two of these (HD 34445 and HD 156668) show strong linear correlation between
their RV trends and S-values, which we can treat as a proxy for magnetic
activity. Therefore, only 2 out of 65 single giant systems show RV evidence
for a stellar binary companion (HD 195019 and HD 145934), and neither of these
systems host giants with $e\geq 0.2$. HD 195019 has a visually resolved wide-
orbit binary companion (Lu et al., 1987), while HD 145934’s status is unclear.

Alternatively, it is possible that most of the single planets in our sample do
in fact have additional sub-Jovian giant companions located outside 10 au. Our
survey may be missing giant planets beyond this distance, and therefore
misidentifying multi giant systems as single giant systems. As seen in Figure
13, the majority of CLS-detected giant planet pairs have a semi-major axis
ratio less than 10, with a peak ratio frequency around 3. This means that our
coldest single giants around 10–20 au could have companions out around 30–60
au, where they would be undetectable if they were small enough and on circular
orbits. Again, this is disfavored by the giant planet fall-off beyond 8 au
reported in Fulton et al. (2021). Additionally, there could be highly-inclined
giant planets in our sample that evade RV detection and enhance the apparent
population of highly eccentric singles.

Figure 13: Observed semi-major axis ratio distribution of two-giant systems.
The curve is a kernel density estimation of the distribution, while the
vertical ticks show the observed set of axis ratios. Horizontal bar shows the
kernel width.

Where do our lonely and highly eccentric giants come from? Out of 50 CLS
single giant planets beyond 0.3 au, 11 have median eccentricities greater than
0.5. CLS contains 30 multi giant systems with at least one giant beyond 0.3
au. If the eccentric singles all used to be scattered multis, that would imply
that 11 out of 41 multi systems, or roughly a quarter, scattered themselves
into singledom. Carrera et al. (2019) found that the observed eccentricity
distribution of giants within 5 au is compatible with scattering accounting
for all giants with $e\geq 0.3$. This supports the claim that many highly
eccentric singles used to be multis that experienced scattering events.
However, our poorly constrained metallicity distribution for high-eccentricity
giant hosts limits our ability to make claims about their original
multiplicity.

We leave this discussion with more questions than when we began, but are
hopeful that future work will elucidate these mysteries. Longer RV baselines
for eccentric single systems and imaging campaigns to look for distant stellar
or substellar companions could clarify the statuses of our single and
eccentric giants, while synthetic population models might have something
interesting to say about the origins of our observed eccentricity and
metallicity distributions.

###  5.2 Ultra-hot Jupiter pile-up

Knutson et al. (2014) performed an RV-trend search for massive companions to
hot Jupiters, and found that $50\pm 10\%$ of giants within 0.7–11 days host an
outer companion between 1–13 $M_{\textrm{J}}$ and 1–20 au. More than two-
thirds of their hot Jupiters with companions have orbital periods within 6
days. In the CLS sample, 3 out of 15 hot Jupiters within 10 days host an outer
companion in the same range (HD 187123, HD 217107, HD 9826). That gives a
lower bound on the companion rate of $20\pm 12\%$.

These two rates are compatible within $2.14\sigma$, and therefore consistent
with the companion occurrence rate reported in Knutson et al. (2014). In order
to make a more direct comparison we would need to calculate the completeness-
corrected companion rate for the hot Jupiters in the Legacy sample, but this
is outside the scope of the current study. It is also worth noting that many
of the companions reported in Knutson et al. (2014) were detected as linear
trends with 5-8 year baselines. These systems would benefit from additional RV
follow-up in order to more rigorously explore potential false positive
scenarios due to stellar activity and better constrain the masses and orbital
semi-major axes of the companions.

###  5.3 Intra-system mass uniformity

Previous studies of sub-Neptune-sized transiting planets have found that
planets located in the same system have correlated masses and radii (Weiss et
al., 2018; Millholland et al., 2017). Our study reveals that giant planets
forming in the same system also appear to have correlated masses. This could
imply that the global properties of protoplanetary disks strongly influence
the final masses of their nascent giant planets, and that these properties are
more important than solid surface density and other properties that vary
radially or over time. That condition could be initial protoplanetary disk
mass, since disk mass determines how much material is available during giant
planet formation, solid mass for large core formation, or disk lifetime, since
quickly dissipating disks have less time to make giants.

Weiss et al. (2018) measured a Pearson correlation of $R=0.65$ for the paired
planet radii of a larger sample than the CLS. This means that the magnitude of
the ‘small peas in a pod’ effect is significantly greater than that of the
‘giant peas in a pod’ effect. The final masses of giant planets are very
sensitive to the relative timing of the onset of runaway accretion. Giant
planets can also cut of the flow of gas to the inner disk, which could affect
the ability of an inner companion to accrete a massive gas envelope. This
might explain the weaker correlation between giant masses.

##  6 Conclusions

We compared and contrasted single giant systems and multi giant systems, and
found several key differences. When we look at giant occurrence with respect
to orbital separation, giant singles have a super-hot Jupiter pile-up not
shared by giant multis, limited by lack of sensitivity to giant planets beyond
20 au. We reproduced previously seen eccentricity distributions for these two
groups, and found that the single giants have a circular orbit pile-up and
long eccentric tail, while the multis are spread across moderate
eccentricities and limited to $e\leq 0.7$. We also found that multi giant
stellar hosts are more metal-rich than single giant hosts, and that total
giant planet mass (inferred from $M\sin i$ ) does not correlate with host
metallicity. Taken together, these results may imply that more metal-rich
stars form more discrete giant planets, possibly due to formation of more
rocky cores. We also found that at least a third of stars that host giant
planets host two or more of them. This implies that giant planet pairs
probably have correlated masses, perhaps that the shared formation environment
of giants influences their final mass limits.

Although this work provides new information about the underlying population of
giants, it also raises unresolved questions about why we see these patterns.
What does the single super-hot Jupiter pile-up mean for theories of hot
Jupiter formation? Are lonely and highly eccentric giants products of planet-
planet scattering? And why are giant planet pair masses correlated? These
questions can be the starting points for new theoretical inquiries into planet
formation.

L.J.R. led the construction of this paper, including performing all analysis,
generating all of the figures, and writing this manuscript. H.A.K. and A.W.H.
advised substantially on the scientific direction of this work. B.J.F.
collaborated on occurrence analysis. A.W.H. acknowledges NSF grant 1753582.
H.A.K. acknowledges NSF grant 1555095. We thank Fei Dai, Lauren Weiss, and
Konstantin Batygin for valuable discussions. 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

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Blunt et al. (2019) Blunt, S., Endl, M., Weiss, L. M., et al. 2019, AJ, 158, 181
  * Bowler et al. (2020) Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63
  * Bryan et al. (2019) Bryan, M. L., Knutson, H. A., Lee, E. J., et al. 2019, AJ, 157, 52
  * Bryan et al. (2016) Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89
  * Carrera et al. (2019) Carrera, D., Raymond, S. N., & Davies, M. B. 2019, A&A, 629, L7
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580
  * Dawson & Johnson (2018) Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175
  * Dawson & Murray-Clay (2013) Dawson, R. I., & Murray-Clay, R. A. 2013, ApJ, 767, L24
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14
  * Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15
  * Huang et al. (2016) Huang, C., Wu, Y., & Triaud, A. H. M. J. 2016, ApJ, 825, 98
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jackson et al. (2021) Jackson, J. M., Dawson, R. I., Shannon, A., & Petrovich, C. 2021, AJ, 161, 200
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 434, L51
  * Knutson et al. (2014) Knutson, H. A., Fulton, B. J., Montet, B. T., et al. 2014, ApJ, 785, 126
  * Lee (2019) Lee, E. J. 2019, ApJ, 878, 36
  * Lissauer et al. (2009) Lissauer, J. J., Hubickyj, O., D’Angelo, G., & Bodenheimer, P. 2009, Icarus, 199, 338
  * Lu et al. (1987) Lu, P. K., Demarque, P., van Altena, W., McAlister, H., & Hartkopf, W. 1987, AJ, 94, 1318
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Millholland et al. (2017) Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441
  * Ogilvie (2014) Ogilvie, G. I. 2014, ARA&A, 52, 171
  * Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273
  * Petrovich & Tremaine (2016) Petrovich, C., & Tremaine, S. 2016, ApJ, 829, 132
  * Rosenthal et al. (2021a) Rosenthal, L. J., Knutson, H. A., Chachan, Y., et al. 2021a, arXiv e-prints, arXiv:2112.03399
  * Rosenthal et al. (2021b) Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021b, ApJS, 255, 8
  * Tange (2011) Tange, O. 2011, ;login: The USENIX Magazine, 36, 42.  http://www.gnu.org/s/parallel
  * Thorngren et al. (2016) Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
  * Van Eylen et al. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48
  * Wright et al. (2009) Wright, J. T., Upadhyay, S., Marcy, G. W., et al. 2009, ApJ, 693, 1084
  * Zhu & Wu (2018) Zhu, W., & Wu, Y. 2018, AJ, 156, 92