# Intrinsic Alignment as an RSD Contaminant in the DESI Survey

Claire Lamman ,1 Daniel Eisenstein,1 Jessica Nicole Aguilar,2 David Brooks,3
Axel de la Macorra,4Peter Doel,3 Andreu Font-Ribera ,5 Satya Gontcho A
Gontcho,2 Klaus Honscheid,6,7 Robert Kehoe,8Theodore Kisner ,2 Anthony Kremin
,2 Martin Landriau ,2 Michael Levi ,2 Ramon Miquel,5,9John Moustakas ,10
Nathalie Palanque-Delabrouille ,2,11 Claire Poppett,2,12,13 Michael
Schubnell,14,15 Gregory Tarlé 15

1Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
2Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720,
USA
3Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK
4Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México
C.P. 04510, México
5Institut de Física de’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain
6Center for Cosmology and AstroParticle Physics, The Ohio State University,
191 West Woodruff Avenue, Columbus, OH 43210, USA
7Department of Physics, The Ohio State University, 191 West Woodruff Avenue,
Columbus, OH 43210, USA
8Department of Physics, Southern Methodist University, 3215 Daniel Avenue,
Dallas, TX 75275, USA
9Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys,
23, 08010 Barcelona, Spain
10Department of Physics and Astronomy, Siena College, 515 Loudon Road,
Loudonville, NY 12211, USA
11IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
12Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way,
Berkeley, CA 94720, USA
13University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA
94720, USA
14Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA
15University of Michigan, Ann Arbor, MI 48109, USA
E-mail<EMAIL_ADDRESS>0000-0002-6731-9329 0000-0002-3033-7312
0000-0003-3510-7134 0000-0001-6356-7424 0000-0003-1838-8528
0000-0003-1887-1018 0000-0002-2733-4559 0000-0003-3188-784X
0000-0003-1704-0781

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We measure the tidal alignment of the major axes of Luminous Red Galaxies
(LRGs) from the Legacy Imaging Survey and use it to infer the artificial
redshift-space distortion signature that will arise from an orientation-
dependent, surface-brightness selection in the Dark Energy Spectroscopic
Instrument (DESI) survey. Using photometric redshifts to down-weight the
shape-density correlations due to weak lensing, we measure the intrinsic tidal
alignment of LRGs. Separately, we estimate the net polarization of LRG
orientations from DESI’s fiber-magnitude target selection to be of order
$10^{-2}$ along the line of sight. Using these measurements and a linear tidal
model, we forecast a 0.5% fractional decrease on the quadrupole of the 2-point
correlation function for projected separations of 40-80 $h^{-1}{\rm Mpc}$. We
also use a halo catalog from the Abacus Summit cosmological simulation suite
to reproduce this false quadrupole.

###### keywords:

methods: data analysis –cosmology: observations – large-scale structure of
Universe – – cosmology: dark energy

††pubyear: 2022††pagerange: Intrinsic Alignment as an RSD Contaminant in the
DESI Survey–B.2 Figure 1: A cartoon demonstrating how an aperture-based
selection can combine with tidal alignment to affect measurements of the
underlying density. Elliptical galaxies will have the maximum concentration of
light on the sky when their primary axis is pointed at the observer. In this
case, more of the light falls within an aperture and it is more likely to be
included in DESI’s fiber magnitude selection. The cartoon on the right shows
galaxies with maximum tidal alignment lying along density filaments which are
parallel to the LOS. These filaments are represented with a blue gradient.
Galaxies in filaments tend to be oriented in the direction of the filament,
and ones between tend to point toward the higher density regions. In this
case, DESI is more likely to select galaxies in denser regions, resulting in
an amplification of this density mode. The opposite effect happens for
filaments which are perpendicular to the LOS (not shown here, see Figure 1 of
Martens et al. (2018)). Since DESI is more likely to select galaxies in
filaments which lie along the LOS, and less likely to select ones in
perpendicular filaments, an anisotropic clustering arises and biases the RSD
signal.

##  1 Introduction

Redshift-Space Distortions (RSD), are an effect often used for measuring the
growth of large-scale structure. On the scale of galaxy clusters, peculiar
velocities of galaxies “smear" structure along the line of sight (LOS) in
redshift space (Jackson, 1972). On larger scales, material falling into over-
dense regions creates a “squashing" effect along the LOS (Kaiser, 1987). The
difference in clustering along versus transverse to the LOS can be described
by the quadrupole of the correlation function, $\xi_{2}$. This needs to be
corrected for to map galaxies in real space, and on large scales is a direct
measurement of the growth rate of structure and can be used to test gravity.

To fully utilize RSD measurements in large spectroscopic galaxy surveys, one
of their important biases must be understood: intrinsic galaxy alignment (IA).
The primary axis of galaxies can be intrinsically aligned with each other (II
correlation) and with the underlying density, or tidal field (GI correlation).
When a galaxy survey has an orientation-dependent selection bias and galaxy
orientations are also correlated with the tidal field, $\xi_{2}$ is directly
affected.

Hirata (2009) used linear models of tidal alignment and orientation-dependent
selections to predict that GI correlations could affect RSD measurements by as
much as 10%. This effect is highly survey-dependent due to its strong
dependence on survey selection and the differences in tidal alignments between
galaxy samples. Martens et al. (2018) and Obuljen et al. (2020) have measured
an anisotropic galaxy assembly bias in the Baryon Oscillation Spectroscopic
Survey (BOSS). Since the velocity dispersion of elliptical galaxies is non-
isotropic and may correlate with axis orientation and tidal environment, this
effect could be a manifestation of the effect described by Hirata (2009). On
the other hand, Singh et al. (2021) followed a similar method to Martens et
al. (2018) and found the Fundamental Plane of BOSS galaxies to be dominated by
systematics and poorly correlated with IA, resulting in a null detection of
the RSD IA bias for BOSS.

As a Stage IV survey, it is necessary to not only detect, but quantify these
biases for the Dark Energy Spectroscopic Instrument (DESI). DESI is in the
midst of a 5-year survey, measuring spectra of over 40 million galaxies within
16,000 deg2 of the sky (DESI Collaboration et al., 2016; Abareshi et al.,
2022).

Successful inference of a galaxy’s spectroscopic redshift depends on target
surface brightness. This is especially true for a large survey like DESI,
which prioritizes survey speed at the cost of higher signal-to-noise. To
impose this explicitly, DESI adopts a surface brightness-dependent cut:
limiting the magnitude within an aperture instead of the objects’ total
magnitude. While this mitigates systematic errors related to surface
brightness, it creates a bias in the 3D orientation of galaxies. Galaxies with
a pole-on orientation have a higher surface brightness and are more likely to
be selected. Since galaxies with tidal alignments tend to point towards
regions of higher density, this can also mean preferentially selecting
galaxies which lie in filaments along the LOS. This results in an enhancement
of clustering in the radial direction and suppression in the transverse
direction, mimicking RSD. The key piece of modeling this effect is relating
the polarizability of the surface brightness selection to the shape of
galaxies viewed from "the side", i.e. transverse to the LOS. This depends on
the details of the light profiles and triaxial shapes of the galaxies (Figure
1).

About 20% of DESI’s targets are Luminous Red Galaxies (LRGs), which fall in
the redshift range $0.4-1.0$ (Zhou et al., 2021). These high-mass, relatively
inactive galaxies exhibit large tidal alignments (Hirata et al., 2007) and are
more affected by an aperture-based selection because they have larger angular
sizes than Emission Line Galaxies (ELGs). Therefore, we chose to focus our
investigation on LRGs as the DESI sample most likely to be substantially
biased by these alignments, although our methods would also work for ELGs.

The two effects that combine to create this bias, GI alignment and selection-
induced polarization, can both be estimated and used to calibrate the
quadrupole $\xi_{2}$. Here we measure the shape-density correlation of LRGs as
projected on the plane of the sky using shapes from the DESI Legacy Imaging
Survey (Dey et al., 2019). We isolate the signal of intrinsic positions from
weak lensing via photometric redshifts, model DESI’s orientation-dependent
selection function, and put our detection in context of $\xi_{2}$ via a linear
tidal model. As an additional test, we use the AbacusSummit cosmological
simulations to reproduce an aperture-based selection and measure the effect on
$\xi_{2}$.

##  2 DESI Catalogs

###  2.1 Imaging

Our measurements of GI alignment were made with LRGs from the Legacy Imaging
Survey, DR9 (Dey et al., 2019). This is the catalog DESI uses to select its
targets, and contains imaging in three bands (g, r, and z) and projected
shapes for sources in $14,000\deg^{2}$ of the extra galactic sky. It also
includes photometry from the Wide-field Infrared Survey Explorer, which
contains $r$ and $W1$ fluxes that are corrected for Milky Way extinction.

The source of each target (after PSF deconvolution) is modeled as several
light profiles at the pixel level using Tractor (Lang et al., 2016). Based on
the fits’ $\chi^{2}$ values, we used shape parameters from the best fit out of
these models: exponential disk, de Vaucouleurs, and Sersic. This is different
from DESI’s default selection, which includes PSF and round-exponential fits,
and a marginalized $\chi^{2}$ criteria to avoid over-fitting bright targets as
round-exponentials. These models were avoided for our measurements, as circles
have no distinguishable orientation.

Quality cuts were applied to target declinations $\delta>-30^{\circ}$ and
galactic latitudes $b>20^{\circ}$. $r-W1$ color correlates well with redshift,
so we used this color for the pair selection and weighting scheme detailed in
3.2. To conform with these weights, color outliers were removed by requiring
$1<r-W1<4.5$. Our final sample contained 17.5 million LRGs.

###  2.2 Spectroscopy

We calibrate our photometric redshifts using a large sample of spectroscopic
redshifts from the DESI Survey Validation (SV) observations. SV is designed to
represent the full survey and is used to assess DESI’s target selection. We
use DESI’s internal SV catalog, Fuji, which comprises of quality observations
taken from 14 December 2020 through 10 June 2021. From this we selected
133,924 LRGs with colors $0.6<r-z$ and $1.5<r-W1<4.5$, and redshifts
$0.001<\text{z}<1.4$.

##  3 Intrinsic Alignment Signal

###  3.1 Alignment Formalism

Figure 2: How angles are defined in our alignment metric. For a given pair of
galaxies, $\epsilon_{1}^{\prime}$ is maximum when the axis of B is parallel to
the separation vector between it and A (most aligned), and minimum when it’s
perpendicular (anti-aligned). $\epsilon^{\prime}_{1}$ also scales with axis
ratio; it approaches 0 as the shape of B becomes more circular.

The projected alignment of galaxies on the sky is quantified with a relative
complex ellipticity (Figure 2). This measures the degree to which a galaxy is
aligned with, and stretched along, a separation vector between it and another
galaxy. Measuring this as a function of the separation vector’s magnitude for
many galaxy pairs is a way to quantify the alignment of LRGs to the underlying
tidal field.

Here, 2D galaxy shapes are modeled as ellipses with a complex ellipticity

$\epsilon=\frac{a-b}{a+b}\exp{2i\phi}$ (1)

where $a$ and $b$ are the primary and secondary axis of the 2D ellipse, and
$\phi$ is the orientation angle of the primary axis, measured East of North.
We define the ellipticity of a galaxy $B$ relative to another galaxy $A$ using
the difference between $B$’s orientation angle, $\phi_{B}$, and its position
angle relative to $A$, $\theta_{BA}$, also measured East of North.

$\theta^{\prime}=\phi_{B}-\theta_{BA}$ (2)

This gives us a relative ellipticity, for which we measure the real component.

$\epsilon^{\prime}=\frac{a_{B}-b_{B}}{a_{B}+b_{B}}\exp{2i\theta^{\prime}}\\\ $
(3)
$\epsilon^{\prime}_{1}=\text{Re}(\epsilon^{\prime})=|\epsilon^{\prime}|\cos{2\theta^{\prime}}\\\
$ (4)

This measurement is averaged over many pairs of galaxies at varying
separations to obtain $w_{\times}(R)$, the 2D shape-density correlation.

Figure 3: Visualization of the matrix used for weighting pairs of galaxies in
the alignment signal based on their color difference. This weighting scheme
was created using a sub-sample of DESI galaxies with measured redshifts and
favors pairs which are more likely to be physically associated with each
other.

###  3.2 Color Weighting

As our signal is a function of transverse separation, the main source of its
dilution is from pairs of galaxies with large radial separations. At the time
of this paper, we do not have spectra for all of the imaged galaxies and so
use color as a redshift proxy. To give pairs which are more likely to be
physically associated a higher weight in the alignment signal, we created a
weighting scheme based off of their $r-W1$ colors.

Using the redshifts DESI has measured so far, described in Section 2.2, we
separated galaxies into 20 bins of $r-W1$ color. We then calculated the
fraction of galaxies which are radially separated by less than $10\ {\rm Mpc}$
for every combination of two colors, based on their redshift difference and
assuming the Hubble flow. The resulting lookup matrix was then used as a
weight when averaging the alignment signal from individual pairs (Figure 3).

###  3.3 Intrinsic Alignment Measurement

The catalog was divided into 10 groups based declination and then each of
those into 10 groups based on right ascension, resulting in 100 sky regions
with an equal number of galaxies in each, 1.8 million. We measured the
projected alignment of neighboring galaxies relative to each galaxy in each
region. This was averaged over 20 bins of angular separation, resulting in 100
determinations of the IA signal. The average and standard error of these 100
measurements at each separation is our projected IA measurement,
$w_{\times}(R)$111code available here: github.com/cmlamman/ellipse_alignment.

Our final determination of $w_{\times}(R)$ for DESI LRGs is displayed in
Figure 4. This signal broadly agrees with our measurement of projected IA in
the Abacus Mock from Section 4, which did not include any misalignments from
the original halo orientations. The similarity between the alignment in LRGs
and raw halo shapes is likely a coincidence due to two opposing effects: halo
orientations are more aligned with the underlying density, which increases
$w_{\times}$, but are rounder than LRGs, which dilutes $w_{\times}$. The LRG
measurements of $w_{\times}$ in each angular bin are statistically independent
of each other, as demonstrated by the covariance of our final $w_{\times}(R)$
signal between the 20 bins of radial separation (Figure 5).

Figure 4: Our final measurement of the projected shape-density correlation of
DESI LRGs (red), which includes weighting based on the color difference in
galaxy pairs. This is compared to the same measurement made with Abacus halos
(dark blue). The light blue line shows the alignment of Abacus halos once the
distribution of their projected shapes was adjusted to match the LRGs, but
does not include galaxy-halo misalignments. Figure 5: The reduced covariance
matrix of $w_{\times}$ between bins of radial separation for our IA
measurement; the identity matrix has been subtracted from this plot. This
demonstrates that the measurements of $w_{\times}$ in each angular bin are
statistically independent of each other.  Figure 6: The shape-density
correlation of DESI LRGs, resulting from measuring the shape of one "neighbor"
galaxy relative to the separation vector between it and another. No color
weights were used for these measurements. Instead, measurements were made
using different color-based restrictions for each pair of galaxies in order to
explore the effects of weak lensing on the IA signal. The dark red line is the
resulting signal when only measuring galaxy pairs which have a very similar
$r-W1$ color, to approximate physical proximity. The orange and yellow lines
are both measurements made on pairs of galaxies which have a large difference
in $r-W1$ color, to emulate pairs which have no physical association. For the
measurement shown in orange, we used pairs in which the neighbor galaxy was
more blue than the other. Therefore we only measured the shape of galaxies
relative to ones behind it, so their shapes were broadly unaffected by weak
lensing. The converse was applied for the signal shown in yellow; here we only
measured the shape of a galaxy if it was much redder than its counterpart.
This means that the measured shape correlation is dominated by weak lensing.

###  3.4 Weak Lensing

Besides intrinsic alignments, the main effect impacting the alignment signal
is gravitational weak lensing. If the shape of a neighbor galaxy is measured
relative to a foreground central galaxy, the neighbor’s light can be
gravitationally distorted by the central before reaching us.

Since we only measure the shape of one galaxy in each pair, weak lensing is
only present when the measured galaxy is in front of the central one.
Therefore a simple way to isolate the weak lensing and IA signals is to set
restrictions on the radial separations of pairs. Using $r-W1$ color again as a
distance proxy, we measured the alignment for sets of pairs with various color
restrictions (Figure 6). The signal from only measuring the shapes of galaxies
relative to their closest color neighbors is comparable to our measurement
using color weighting. The signal from only measuring galaxy shapes relative
to background galaxies is consistent with $0$ above separations of a few Mpc,
and the signal from measuring galaxy shapes relative to foreground galaxies
is, as expected, opposite in sign to the intrinsic signal.

To check whether the lensing signal is consistent with expectations, we
consider the following approximate model. The net effect of weak lensing acts
in opposition to the IA signal, as it creates a tangential shear on the sky:

$\gamma_{t}=\frac{\bar{\Sigma}(<r)-\Sigma(r)}{\Sigma_{\rm crit}}$ (5)

where $\bar{\Sigma}(<r)$ is the average surface density with some transverse
distance $r$. Here, $r_{0}=7.78\ {\rm Mpc}/h$ is the 3D correlation length for
DESI clustering (Kitanidis et al., 2020), $\beta=2.15$ is the clustering bias
for DESI LRGs (Zhou et al., 2021), and $\rho_{0}=2.68\times 10^{-30}\text{ g
}\text{cm}^{-3}$ is the critical matter density of the Universe from Planck
2018 (Collaboration et al., 2020). $D_{S}$, $D_{L}$, and $D_{L}S$ are the
distances to the source, distance to the lens, and distance between them,
respectively. Assuming the correlation function goes as
$\xi(r)=r_{0}^{2}\rho_{0}/\beta r^{2}$,

$\bar{\Sigma}(<r)=\frac{2\pi}{r}r_{0}^{2}\frac{\rho_{0}}{\beta}$ (6)

$\Sigma(r)$ is the average surface density at $r$

$\Sigma(r)=\frac{r_{0}^{2}}{r}\pi\frac{\rho_{0}}{\beta}$ (7)

and $\Sigma_{\rm crit}$ is the critical mean density, above which the light of
a source is split into multiple images.

$\Sigma_{crit}=\frac{c^{2}}{4\pi
G}\frac{D_{S}}{(1+z_{l})D_{L}D_{LS}}\approx\frac{c^{2}D_{S}}{4\pi
GD_{L}D_{LS}}$ (8)

To connect this to our alignment formalism described in Section 3.1, the
tangential shear is defined as

$\gamma_{t}=\frac{a-b}{a+b}e^{2i\phi}$ (9)

where $\phi$ is the azimuthal angle of the source galaxy’s primary axis with
respect to the lens. This results in the relation

$\bar{\epsilon_{1}}^{\prime}=\frac{\bar{\gamma}_{t}}{-2}$ (10)

To measure this in our sample, we used photometric redshifts to estimate
$D_{S}/D_{L}D_{LS}$ for every pair of galaxies, and average the result. We
used a simple, linear fit of our DESI spectroscopic sample for redshifts:

$\text{z}=0.25(r-W1)-0.02$ (11)

The resulting lensing estimation is shown in Figure 6 and agrees well with the
IA measurement made when limiting to pairs we expect are only affected by
lensing. This final IA signal is likely still diluted by weak lensing. However
we did not develop a more sophisticated adjustment for lensing, as DESI’s
first year of spectra will allow us to sufficiently isolate physically-
associated pairs.

##  4 IA with Abacus Mock Catalog

To contextualize the measured IA signal, we built a mock catalog from the
AbacusSummit CompaSO halo catalog (Hadzhiyska et al., 2021). AbacusSummit is a
suite of large, high-accuracy, high-resolution cosmological simulations made
with the Abacus N-body code (Maksimova et al., 2021). We used halos from a box
with comoving $2000\ h^{-1}{\rm Mpc}$ sides, simulated at $z=0.725$.

We mapped the halos’ comoving positions to redshift and sky coordinates by
placing an observer $1700\ h^{-1}{\rm Mpc}$ away from the center of the box
along the $x$-axis. To have an even sky distribution and consistent redshift
range across the sky, we set boundaries of $\pm 12^{\circ}$ in right ascension
and declination, with a redshift range of $0.51<\text{z}<0.97$.

We then selected the largest halos to match both the LRG density of our DESI
sample within this redshift range, $7.3\times 10^{-4}\ (h^{-1}{\rm
Mpc})^{-3}$, and the redshift distribution from DESI spectra. Our final mock
catalog contains 766,341 halos.

To imitate the DESI Legacy Imaging Survey colors, we used a catalog of DESI
LRG spectroscopic redshifts, sorted by redshift. For each halo, we identified
the LRG with the closest redshift. We then smoothed our selection by sampling
a neighboring LRG from a Gaussian of indices centered at the index of closest
redshift and with a width of 300. TAfter taking the $r-W1$ color from the LRG
with the resulting index, we again smoothed by sampling a Gaussian centered at
that color, with a width of $\sigma=0.03$. These smoothing parameters
sufficiently reproduced the observed data spread, and variations of them do
not significantly affect the measured alignment signal.

The Abacus 3D halo shapes are modeled as triaxial ellipsoids. A common method
for finding the projected axis ratios of ellipsoids can be found in Binney
(1985). For measuring the alignment of galaxy shapes, we additionally need the
orientation angle of the projected shape. Therefore, we adapted the method
derived in Gendzwill & Stauffer (1981) to project triaxial ellipsoids onto the
celestial sphere. Our process for obtaining the axis ratio and orientation of
an ellipsoid projected along an arbitrary viewing angle can be found in
Appendix A.

Since halos are rounder than LRGs, so we mapped the axis ratios of the
projected halos to the LRG axis ratio distribution. We adjust each axis ratio,
$b/a=d$, with the empirical function:

$d^{\prime}=1+1.1(d-1)-2.035(d-1)^{2}+1.76(d-1)^{3}$ (12)

This function correctly reproduces the number of observed axis ratios which
fall in 100 bins between 0 and 1. We made no adjustments for the orientations
of halos.

Using the same color-weighting scheme as described in Section 3, we measure
the projected shape-density correlation of our resulting mock catalog. The
result can be seen in Figure 4. The higher amplitude is likely due to the
simulation not including the effects of weak lensing and the higher degree of
alignment in halos compared to galaxies. Tenneti et al. (2014) estimates
large, central galaxies at DESI redshifts to be misaligned with their host
halo by an average of around 10-20∘. This propagates to a $w_{\times}$ signal
that is 75-94% of the same signal without misalignments.

##  5 Modeling alignment - $\xi_{2}$ Correlation

###  5.1 Linear Tidal Model

We adopt a linear tidal model to connect IA and DESI’s shape selection bias
with the quadrupole of the correlation function, $\xi_{2}$. This approximation
assumes that the projected shapes of galaxies are linearly related to the
projected density distribution and holds for LRGs above projected separations
of $10\ h^{-1}{\rm Mpc}$ (Catelan & Porciani, 2001; Hirata & Seljak, 2004;
Singh et al., 2015; Troxel & Ishak, 2015).

We define $\nabla^{2}\phi=\rho$, where $\rho$ is the fractional over density.
The tidal tensor is then the traceless combination
$T_{ij}=\partial_{i}\partial_{j}\phi-(1/3)\delta^{K}_{ij}\nabla^{2}\phi$,
where $\delta^{K}$ is the Kronecker delta.

We model the mean 3D ellipticity of a triaxial galaxy as $\tau T_{ij}$, where
the axis lengths behave as $I+\tau T$. For this derivation, we assume that 2D
projections of such galaxies behave as the 2D projection of these lengths. The
mean eccentricity tensor must also be traceless, so for a projection with
$\alpha$, $\beta$ = $\\{x,y\\}$, the projected ellipticity is given as
$\epsilon_{\alpha\beta}=\tau(T_{\alpha\beta}+T_{zz})$, where we used
$T_{xx}+T_{yy}=-T_{zz}$.

Using Fourier-space conventions and the matter power spectrum
$\left<\tilde{\rho}(\vec{q})\tilde{\rho}^{*}(\vec{k})\right>=(2\pi)^{3}P(\vec{k})\delta^{D}(\vec{q}-\vec{k})$,
we have the tidal tensor model:

$\begin{split}T_{ij}(\vec{r})=\left(\partial_{i}\partial_{j}-\frac{\delta^{K}_{ij}}{3}\nabla^{2}\right)\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{\phi}(\vec{k})e^{i\vec{k}\cdot\vec{r}}\\\
=\int\frac{d^{3}k}{(2\pi)^{3}}\left(\frac{k_{i}k_{j}-\delta^{K}_{ij}k^{2}/3}{k^{2}}\right)\tilde{\rho}(\vec{k})e^{i\vec{k}\cdot\vec{r}}\end{split}$
(13)

###  5.2 Shape-Density Correlation

To consider the correlation of galaxy ellipticity with surface density, we
begin with an expression for the projected fractional overdensity for a survey
of functional depth $L$ and uniform mean density $\rho$:

$\Sigma(\vec{R})=\frac{1}{L}\int dz\;\rho(\vec{R},z)$ (14)

where $\hat{z}$ is along the LOS and $\vec{R}$ is projected separation. $L$ is
a measure of how far along the LOS we average when measuring
$\epsilon_{\text{LRG}}$. As our survey is not homogeneous, we generalize $L$
to an expression of $N(z)$. Using the weights we give each galaxy pair $w$, we
sum over all combinations of color bins $B1$, $B2$, and galaxy pairs $j$,
$k$.This is averaged per-bin and multiplied by the depth of that bin $B_{d}$.

$L=B_{\text{d}}\frac{\Sigma_{B1}\Sigma_{B2}\Sigma_{\text{i}}\Sigma_{\text{j}}w(i,j)}{\Sigma_{B1}\Sigma_{\text{i}}\Sigma_{\text{i}}w(i,j)}$
(15)

The projected ellipticity is
$\hat{R}_{\alpha}\epsilon_{\alpha\beta}\hat{R}_{\beta}$. For the average, we
can just consider the $\hat{R}=\hat{x}$ direction. The shape-density
correlation projected onto the plane of the sky is then given as

$w_{\times}(R)=\left<\epsilon_{xx}\Sigma(R\hat{x})\right>=-\frac{\tau}{L}\left<T_{yy}\int
dz\;\rho(R\hat{x},z)\right>$ (16)

We compute and simplify this expression (Appendix B.1) to

$w_{\times}(R)=\frac{\tau}{3L}R^{2}\frac{d}{dR}\left[\frac{1}{R^{2}}\Psi(R)\right]$
(17)

where we introduce

$\Psi(R)=-\frac{d\Phi}{dR}=\int\frac{K\;dK}{2\pi}\frac{P(K)}{K}J_{1}(KR)$ (18)

$K$ is 2D Fourier Space, $P$ is the power spectrum, and $J_{1}$ is the first
Bessel function.

$\tau$ can be inferred from our measurement of the shape-density correlation,
$\left<\epsilon_{xx}\Sigma(R\hat{x})\right>$, showing that the LOS shape and
$\xi_{2}$ are correlated. We estimate $\tau$ as

$\tau_{\rm
obs}=\frac{3Lw_{\times}}{R^{2}\frac{d}{dR}\left[\frac{1}{R^{2}}\Psi\right]},$
(19)

with our IA measurement, $w_{\times}$, and averaging over angular scales $R$.

###  5.3 Shape - $\xi_{2}$ Correlation

Next, we turn to the correlation of shapes with the LOS $\xi_{2}$. To obtain
the relation between 3D shapes and the LOS, we consider shapes viewed from a
direction transverse to the LOS, i.e. an axis perpendicular to the projection
axis above.

Since we define $\xi_{2}$ as
$\xi(r,\mu)=\sum_{\ell}\xi_{\ell}(r)L_{\ell}(\mu)$, with $\mu$ the cosine of
the angle to the LOS, $\xi_{2}$ is given as

$Q(r)=5\int\frac{d^{2}\hat{r}}{4\pi}\rho(\vec{r})L_{2}(\mu)$ (20)

where $\int d^{2}\hat{r}$ indicates the 2D integral over the unit vector
$\hat{r}$.

To express the transversely viewed shape, we take the average of
$\epsilon_{xz}$ and $\epsilon_{yz}$, each after the correction to a traceless
quantity. For the $x$–$z$ projection, we have
$\epsilon_{\alpha\beta}=\tau(T_{\alpha\beta}+T_{yy})$, where the relevant
quantity is $\epsilon_{zz}$. Averaging with the $y$–$z$ projection, we have a
transverse averaged eccentricity

$\epsilon_{zz}=\tau\left(T_{zz}+\frac{T_{xx}+T_{yy}}{2}\right)=\frac{\tau}{2}T_{zz}$
(21)

Considering projections along $\hat{x}\pm\hat{y}$ also yield $T_{zz}/2$ as the
only $m=0$ support.

The details of computing and simplifying $\left<\epsilon_{zz}Q(r)\right>$ can
be found in Appendix B.2, which result in

$\left<\epsilon_{zz}Q(r)\right>=-\frac{\tau}{3}\int\frac{q^{2}dq}{2\pi^{2}}P(q)j_{2}(qr).$
(22)

This expression is averaged over radial bins of the correlation function,
resulting in averages of $j_{2}(qr)$.

###  5.4 Bias on $\xi_{2}$

We expect the mean shape to be elongated along the LOS due to DESI’s target
selection, i.e. a non-zero mean $\epsilon_{zz}$ (Section 6). We call this LOS
polarization $\epsilon_{\rm LRG}$.

Assuming $\epsilon_{zz}$ and the quadrupole signature $Q$ are Gaussian
distributed, correlated, random variables, a non-zero
$\left<\epsilon_{zz}\right>$ will result in a non-zero mean $Q(r)$ as

$\left<Q\right>=\left<\epsilon_{zz}\right>\frac{\left<\epsilon_{zz}Q\right>}{\left<\epsilon_{zz}^{2}\right>}.$
(23)

where the expectation values come from summing over each galaxy. From our
tidal model,

$\displaystyle\left<\epsilon_{zz}^{2}\right>$ $\displaystyle=$
$\displaystyle\frac{\tau^{2}}{4}\left<T_{zz}^{2}\right>$ (24) $\displaystyle=$
$\displaystyle\frac{\tau^{2}}{45}\int\frac{q^{2}dq}{2\pi^{2}}P(q).$ (25)

This is the variance in the density field $\sigma^{2}$, hence
$\left<\epsilon_{zz}^{2}\right>=\tau^{2}\sigma^{2}/45$.

Combining the above results, we obtain an expression for the quadrupole
signature arising from GI alignment and a shape-dependent selection bias:

$\xi_{\rm GI}=\left<Q(r)\right>=\epsilon_{\rm
LRG}\frac{\tau}{3\left<\epsilon_{zz}^{2}\right>}\int\frac{q^{2}dq}{2\pi^{2}}P(q)j_{2}(qr)$
(26)

##  6 Modeling DESI’s selection effects

In Section 3 we measured how the shapes of galaxies projected onto the sky are
aligned with the density field. To infer how this affects RSD measurements, we
need to estimate the extent of DESI’s orientation-dependent selection bias.
Since pole-on galaxies have a higher surface brightness and are more likely to
pass selection, we expect a net orientation of galaxies along the LOS, or
polarization $\epsilon_{LRG}$. This is defined as the ellipticity (Equation 4)
relative to the LOS.

We estimate this by using a parent catalog of LRGs which is similar to the
sample described in Section 2.1, except without the fiber magnitude cut. We
assign each parent LRG a 3D galaxy light profile, then simulate images of each
profile from all viewing angles, without any extinction from internal dust.
The polarization is the average $\epsilon_{\text{1,LOS}}$ of all 3D profiles
which pass selection.

###  6.1 Parent Sample

We estimate polarization using a subsection of DESI LRGs in an area of the sky
with the best-resolved shapes, with right ascension and declination limits of
0h0m0s $<\alpha<$ 0h40m0s and 0${}^{\circ}<\delta<$ 5∘. This is in the South
Galactic Cap (SGC) and part of the Legacy Imaging Survey’s DES region. This
parent sample of 41120 objects has the same criteria as DESI’s final target
selection, except without the fiber z-magnitude cut of
$z^{\prime}_{\text{fiber}}<21.61$ for the SGC. The fiber magnitude comes from
the light within a 1.″5 aperture after convolving the shape model with a
standardized PSF. This somewhat isolates the fiber magnitude from seeing
variations, so we can safely use shapes from an area with the best seeing
without impacting the distribution of underlying shapes. As in Section 2, we
also use shape parameters from the best-fit, non-circular, model.

$0.95\%$ of this sample have the same fiber $z$-magnitude as total
$z$-magnitude. This indicates that these objects are either stars or
unresolved galaxies. We ignored these objects for our analysis, but a more
thorough simulation would involve simulating galaxies through the Tractor
pipeline, as is done with the Obiwan project (Kong et al., 2020).

###  6.2 Light Profiles

Our light profile for each galaxy begin as a realization of 100,000 points.
This representation allows us to rapidly apply the triaxial axis lengths,
rotations, and projections, as well as to apply a 2D Gaussian PSF (by
deflections) and the eventual fiber aperture cut.

The points for a given galaxy are distributed in 3D based on its best-fit
shape model from the parent catalog. DESI’s Tractor pipeline represents
projected galaxy shapes as a mixture of Gaussians (Hogg & Lang, 2013). To de-
project these into 3D profiles, we take advantage of the fact that a 3D
Gaussian projects to a 2D Gaussian. Therefore the 2D Gaussian mixture fits
allow us to immediately construct a 3D model. This was done for all parent
LRGs with a best-fit profile of de Vaucouleurs, exponential, and round-
exponential LRGs. Relatively few LRGs were fit best with a Sersic profile.
These tend to be bright enough that they are not near the aperture magnitude
cut and therefore less affected by this biased selection; for simplicity we
modeled these with a Hernquist profile (Hernquist, 1990).

Figure 7: A comparison of the axis ratios of LRGs in our parent sample and the
projected axis ratios from a distribution of triaxial shapes. These are the
triaxial shapes used in our polarization estimate. The spike at $b/a=1$ in the
parent sample is artificial, likely due to poor shape fitting. Figure 8:
Results from our N-body reproduction of DESI’s target selection. There are two
main flux cuts on the LRGs: a sliding $r-W1$ vs $W1$ cut which dominates at
bluer colors, and the $z_{\text{fiber}}$ cut which dominates at redder colors.
The full parent sample is shown on the left and a closer look near the fiber
magnitude cut on the right. Each galaxy was assigned a triaxial shape, which
was rotated to 100 random orientations. Its polarization is the average
ellipticity relative to the light of sight of the objects which passed an
aperture-magnitude cut. For target selection, we find that the orientations of
shapes matters only for objects very close to the fiber magnitude cut, and is
more likely to matter for more elliptical galaxies.

(a)

(b)

Figure 9: Properties of selected targets from our simulated images. a: the $z$
mag and color taken from a parent sample of LRGs. Each square is colored by
the average polarization in that bin. We attribute the drop in polarization
near the highest $z$ mag to a selection effect: in order for these targets to
pass selection they must have a fiber magnitude very close to their total
magnitude, resulting in more compact shapes and a dampened polarization. b:
the polarization of selected targets binned by color. We expect this trend,
since fainter galaxies tend to fall closer to the fiber-magnitude cut (Figure
8). A higher polarization for redder colors could lead to an increased
$\xi_{2}$ bias at higher redshifts and mimic structure growth.

###  6.3 Polarization Estimate

For each object in the parent sample, we assign a triaxial shape based on its
projected shape. These were randomly drawn from the expected distribution of
triaxial shapes for bright ($M_{r}>-19$), medium ($2<r_{gal}<7h^{-1}kpc$)
ellipticals in imaging from the Sloan Digital Sky Survey (Padilla & Strauss,
2008). 41120 3D shapes were projected along a random viewing angle and ranked
by the axis ratio of the resulting ellipse. The LRG parent sample was also
sorted by axis ratio, and matched with the triaxial shape corresponding to the
projected shape of the same rank.

To test these triaxial shapes, we viewed them each from a different angle and
compared the projected axis ratios to our parent sample (Figure 7). These
distributions are not identical; note the artificial spike in the parent
sample at $b/a=1$ which is likely from poor shape fitting. Differences in the
distributions could also be due to shape-dependent fitting biases in Tractor,
or imperfect distributions from Padilla & Strauss (2008), including shape
evolution from z $=0$ or internal obscuration.

The point positions from Section 6.2 were scaled by the assigned three axis
lengths for each galaxy. They were then rotated to 100 random orientations and
projected along one axis. The resulting ‘images’ were scaled using the ratio
of the observed half-light radius and the average half-light radius of all
model orientations. We next need to emulate an observation in 1″seeing.
Instead of convolving with a Gaussian, we took the quicker approach of adding
pre-computed, 2D deflections to the projected points. The fiber magnitude was
estimated by from the fraction of points which fell within an 1.″5 - diameter
aperture, and the observed total magnitude of the LRG. The light profiles used
did not perfectly replicate the observed $z_{\text{fiber}}$ values, so we
added a calibration factor to the N-body fiber magnitude for each of the four
light profiles to match the true $z_{\text{fiber}}$ median. Objects with a
fiber magnitude less than 21.61 passed selection.

For each simulated image which passed selection, we measured the corresponding
3D profile’s complex ellipticity relative to the LOS. This is the same
convention as Equation 4, except shapes are projected in the transverse
direction. The average of these is our polarization $\epsilon_{\rm LRG}$.
54.2% of our simulated galaxy images passed the fiber magnitude cut, similar
to the actual value of 52.9%. The polarization for these galaxies is
$0.0087\pm 0.0002$. By determining the selection of a set of orientations for
each galaxy shape, we can also estimate which galaxies in the original sample
may have an orientation-dependent selection (Figure 8). To see what
polarization DESI can expect in its targets, we’ve plotted the average
polarization in bins of $z$ mag and $r-W1$ color (Figure 9a).

We also find that the redder LRGs may be more affected by orientation. This
translates to a correlation between redshift and polarization, which could
affect studies of structure evolution (Figure 9b).

Figure 10: Estimated impact on DESI’s measurement of the RSD quadrupole. This
about corresponds to a 0.5% error for separations around 40-80 $h^{-1}{\rm
Mpc}$ based on SDSS-III BOSS measurements (Anderson et al., 2014). The
precision of our estimate is shown in the shaded region as the standard error.

##  7 Estimate of $\xi_{2}$ Bias

At this point, we have measured all the necessary components to estimate the
$\xi_{2}$ signature arising from IA and DESI’s selection bias.

$\epsilon_{\rm LRG}$, the polarization of galaxy shapes along the LOS, is
measured in Section 6. $\langle\epsilon_{zz}^{2}\rangle$ is the variance of
the real part of the complex ellipticities which describe the shapes of DESI’s
LRGs and is 0.031. We used the power spectrum, $P(K)$ from AbacusSummit
(Maksimova et al., 2021).

$\tau$ is a function of effective depth $L$, or how far along the line of
sight we average when measuring $\epsilon_{\text{L}RG}$. This was estimated
using the color weighting scheme from Section 3.2 to be $L=865h^{-1}{\rm Mpc}$
It also depends on the projected shape-density correlation of LRGs $w_{\rm
x}(r)$ which we measured in Section 3. Averaging over the bins of projected
separation, we estimate $\tau_{\rm obs}=-0.174$.

Using equations 26 and 19 to bring everything together, we determine
$r^{2}\xi_{\rm GI}$ to be $3.5h^{-1}$Mpc around 10-80 $h^{-1}{\rm Mpc}$. The
full separation dependence is shown in Figure 10. SDSS-III measures
$r^{2}\xi_{2}$ at these scales to be near 75 ($h^{-1}$Mpc)2 (Anderson et al.,
2014). This puts our estimate of the fractional error on $\xi_{2}$ around 0.5%
at 40-80 $h^{-1}$Mpc.

##  8 $\xi_{2}$ Bias with Abacus

To demonstrate that an aperture selection produces a false $\xi_{2}$ signature
and test our linear tidal model connecting the GI and RSD signals, we next
model the problem using Abacus Summit simulations.

As in Section 4, we started with a 2000 $h^{-1}$Mpc box of large halos and
mapped their positions to redshift, right ascention, and declination. Sky cuts
were applied to ensure a uniform sky distribution at each redshift. 3D Sersic
profiles of 100,000 points were generated for each halo, as in Sections
6.2-6.3, except using the halo’s original triaxial shape. The half-light
radius used for each halo was drawn from a distribution matching the physical
radii of the DESI LRG parent sample and scaled using the average half-light
radii of the point profile projected to 10 random orientations. We counted the
number of points which fell within a 1.″5 aperture and measured the shape of
each halo projected both on the sky and relative to the LOS.

To see how an aperture selection impacts the $\xi_{2}$ measurement, we created
two samples: one without any selection, and one only with halos containing
more than than 48,000 points within the aperture, which corresponds to 50% of
the halos. We measured $\xi_{2}(r)$ for both sets in real space space and in
redshift space, using the halo’s original velocities.

$\xi_{2}(r)$ was determined using the Landy-Szalay Estimator (Landy & Szalay,
1993) and averaged over 10 sets of randoms, generated with random right
ascension and declinations for each redshift. This entire process was done for
5 Abacus Summit simulation boxes, and their average $\xi_{2}(r)$ and standard
error is shown in Figure 11.

Figure 11: $\xi_{2}$ measurements of Abacus halos catalogs with and without an
aperture-based selection. The top two lines include Redshift-Space
Distortions. The aperture-based selection creates an artificial, non-zero RSD
signature, which acts in opposition to $\xi_{2}$ on large scales.

As in Section 7, we used our linear tidal model to predict the $\xi_{2}$ bias
caused by the aperture selection for this halo catalog. We measured the
projected intrinsic alignment of the halo catalog in radial bins which
resulted in an effective survey depth, $L$, of 138 $h^{-1}$Mpc. The
polarization due to aperture cut was $\epsilon_{LRG}=7.6\pm 0.1\times
10^{-3}$. The resulting prediction is compared to the model in Figure 12.

We expect the bulk of the disagreement between these two simple models to be
due to the linear approximation, which does not hold at lower separations, and
simplifications in the demonstration mock. The largest simplification here is
that every galaxy is modeled with a Hernquist light profile. Any profiles
which are denser than reality will underestimate the polarization due to
aperture selection. However, the Abacus approximation is comparable to the
prediction from the linear model and serves as a adequate demonstration of how
a false $\xi_{2}$ signature can arise for DESI.

Figure 12: The artificial RSD signature induced by an aperture-based
selection. Here we compare the prediction made with our linear tidal model to
the measured difference in $\xi_{2}$ between a halo catalog with and without
the aperture selection. This is the difference between the measurements of
$\xi_{r}$ in Figure 11 without RSD. The model gets close, but appears to be
about 40% lower.

##  9 Conclusion

The objective of this study is to determine the approximate impact on DESI’s
RSD measurements due to an orientation bias in LRGs. We have demonstrated that
the effect is significant for DESI and estimate a 0.5% fractional decrease of
$\xi_{2}$ for separations of 40-80 $h^{-1}{\rm Mpc}$. DESI forecasts a total
$f\sigma_{8}$ around 0.4-0.7% (with ELG and LRGs combined), so it is important
to mitigate this effect.

To reduce the effects of intrinsic alignment for DESI, simple yet severe
choices involve only measuring $\xi_{2}$ in galaxy subsamples, perhaps cut by
total magnitude or color. More practically, our estimate could be used for
calibration.

As the DESI survey progresses and the precision in $\xi_{2}$ increases, there
are several opportunities to improve our bias estimate. Our estimate is
directly proportional to the measured polarization $\epsilon_{\text{LRG}}$ and
IA signal $w_{\text{x}}$, both of which include systematic uncertainties. The
main systematic uncertainty in our polarization estimate arises from the
choice of the triaxial shape distribution. We expect the majority of our
galaxies to be prolate (Padilla & Strauss, 2008), which are more affected by
selection bias than oblate and result in a higher polarization. We match the
expected distribution of projected shapes in a region of the sky with the best
shape fits, but 5.6% of galaxies in this subsample are fit as circles,
creating an artificial spike at $b/a=1$ (Figure 7). A better estimate could be
made with more accurate shapes, ie from the Dark Energy Survey (Gatti et al.,
2021), or reproducing Padilla & Strauss (2008) with DESI’s LRGs.

Although partially mitigated by color weighting, the IA signal in this work is
reduced by weak lensing and diluted by the inclusion of pairs which have large
radial separations. This will be drastically improved with DESI’s first year
of data, which contains 2.5 million quality LRG spectra. The LOS distance we
average over due to uncertainty in radial distances, $L=865h^{-1}{\rm Mpc}$,
will decrease by a factor of at least 20 with redshifts. Advancing our ability
to measure IA for only pairs of galaxies which are physically associated will
be the strongest improvement to the false $\xi_{2}$ estimate.

## Acknowledgements

CL wishes to acknowledge Sean Moss for his support, software and otherwise,
and Tanveer Karim for many helpful conversations. The authors also thank Eric
Jullo, Eusebio Sanchez, Robert Kehoe, and Zachary Slepian for their input
through the collaboration review process.

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE1745303, the U.S. Department
of Energy under grant DE-SC0013718, and NASA under ROSES grant
12-EUCLID12-0004. DJE is further supported by the Simons Foundation.

This research is also supported by the Director, Office of Science, Office of
High Energy Physics of the U.S. Department of Energy under Contract No.
DE–AC02–05CH11231, and by the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility under the same contract;
additional support for DESI is provided by the U.S. National Science
Foundation, Division of Astronomical Sciences under Contract No. AST-0950945
to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the
Science and Technologies Facilities Council of the United Kingdom; the Gordon
and Betty Moore Foundation; the Heising-Simons Foundation; the French
Alternative Energies and Atomic Energy Commission (CEA); the National Council
of Science and Technology of Mexico (CONACYT); the Ministry of Science and
Innovation of Spain (MICINN), and by the DESI Member Institutions:
https://www.desi.lbl.gov/collaborating-institutions.

The DESI Legacy Imaging Surveys consist of three individual and complementary
projects: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona
Sky Survey (BASS), and the Mayall z-band Legacy Survey (MzLS). DECaLS, BASS
and MzLS together include data obtained, respectively, at the Blanco
telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok
telescope, Steward Observatory, University of Arizona; and the Mayall
telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the
Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation. Pipeline
processing and analyses of the data were supported by NOIRLab and the Lawrence
Berkeley National Laboratory. Legacy Surveys also uses data products from the
Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), a project of
the Jet Propulsion Laboratory/California Institute of Technology, funded by
the National Aeronautics and Space Administration. Legacy Surveys was
supported by: the Director, Office of Science, Office of High Energy Physics
of the U.S. Department of Energy; the National Energy Research Scientific
Computing Center, a DOE Office of Science User Facility; the U.S. National
Science Foundation, Division of Astronomical Sciences; the National
Astronomical Observatories of China, the Chinese Academy of Sciences and the
Chinese National Natural Science Foundation. LBNL is managed by the Regents of
the University of California under contract to the U.S. Department of Energy.
The complete acknowledgments can be found at https://www.legacysurvey.org/.

The authors are honored to be permitted to conduct scientific research on
Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the
Tohono O’odham Nation."

## Data Availability

The DESI Legacy Imaging Survey is publicly available at legacysurvey.org.
AbacusSummit simulations are publicly available at abacusnbody.org. Code for
projecting ellipsoids and generating light profiles can be found at
github.com/cmlamman/ellipse_alignment.

All data plotted in this paper are available at zenodo.org/record/7058448

## References

  * Abareshi et al. (2022) Abareshi B., et al., 2022, Technical report, Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument, https://ui.adsabs.harvard.edu/abs/2022arXiv220510939A. Lawrence Berkeley National Laboratory, https://ui.adsabs.harvard.edu/abs/2022arXiv220510939A
  * Anderson et al. (2014) Anderson L., et al., 2014, Monthly Notices of the Royal Astronomical Society, 441, 24
  * Binney (1985) Binney J., 1985, Monthly Notices of the Royal Astronomical Society, 212, 767
  * Catelan & Porciani (2001) Catelan P., Porciani C., 2001, Monthly Notices of the Royal Astronomical Society, 323, 713
  * Collaboration et al. (2020) Collaboration P., et al., 2020, Astronomy and Astrophysics, 641, A6
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, Technical report, The DESI Experiment Part II: Instrument Design, https://ui.adsabs.harvard.edu/abs/2016arXiv161100037D. Lawrence Berkeley National Laboratory, https://ui.adsabs.harvard.edu/abs/2016arXiv161100037D
  * Dey et al. (2019) Dey A., et al., 2019, The Astronomical Journal, 157, 168
  * Gatti et al. (2021) Gatti M., et al., 2021, Monthly Notices of the Royal Astronomical Society, 504, 4312
  * Gendzwill & Stauffer (1981) Gendzwill D. J., Stauffer M. R., 1981, Journal of the International Association for Mathematical Geology, 13, 135
  * Hadzhiyska et al. (2021) Hadzhiyska B., Eisenstein D., Bose S., Garrison L. H., Maksimova N., 2021, Monthly Notices of the Royal Astronomical Society, 509, 501
  * Hernquist (1990) Hernquist L., 1990, The Astrophysical Journal, 356, 359
  * Hirata (2009) Hirata C. M., 2009, Monthly Notices of the Royal Astronomical Society, 399, 1074
  * Hirata & Seljak (2004) Hirata C. M., Seljak U., 2004, Physical Review D, 70, 063526
  * Hirata et al. (2007) Hirata C. M., Mandelbaum R., Ishak M., Seljak U., Nichol R., Pimbblet K. A., Ross N. P., Wake D., 2007, Monthly Notices of the Royal Astronomical Society, 381, 1197
  * Hogg & Lang (2013) Hogg D. W., Lang D., 2013, Publications of the Astronomical Society of the Pacific, 125, 719
  * Jackson (1972) Jackson J. C., 1972, Monthly Notices of the Royal Astronomical Society, 156, 1P
  * Kaiser (1987) Kaiser N., 1987, Monthly Notices of the Royal Astronomical Society, 227, 1
  * Kitanidis et al. (2020) Kitanidis E., et al., 2020, Monthly Notices of the Royal Astronomical Society, 496, 2262
  * Kong et al. (2020) Kong H., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3943
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, The Astrophysical Journal, 412, 64
  * Lang et al. (2016) Lang D., Hogg D. W., Mykytyn D., 2016, Astrophysics Source Code Library, p. ascl:1604.008
  * Maksimova et al. (2021) Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, Monthly Notices of the Royal Astronomical Society, 508, 4017
  * Martens et al. (2018) Martens D., Hirata C. M., Ross A. J., Fang X., 2018, Monthly Notices of the Royal Astronomical Society, 478, 711
  * Obuljen et al. (2020) Obuljen A., Percival W. J., Dalal N., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 058
  * Padilla & Strauss (2008) Padilla N. D., Strauss M. A., 2008, Monthly Notices of the Royal Astronomical Society, 388, 1321
  * Singh et al. (2015) Singh S., Mandelbaum R., More S., 2015, Monthly Notices of the Royal Astronomical Society, 450, 2195
  * Singh et al. (2021) Singh S., Yu B., Seljak U., 2021, Monthly Notices of the Royal Astronomical Society, 501, 4167
  * Tenneti et al. (2014) Tenneti A., Mandelbaum R., Di Matteo T., Feng Y., Khandai N., 2014, Monthly Notices of the Royal Astronomical Society, 441, 470
  * Troxel & Ishak (2015) Troxel M. A., Ishak M., 2015, Physics Reports, 558, 1
  * Zhou et al. (2021) Zhou R., et al., 2021, Monthly Notices of the Royal Astronomical Society, 501, 3309

##  Appendix A Projection of Triaxial Ellipsoids

This section details how we obtained the axis ratios and orientations of
projected triaxial ellipsoids for our mock catalogs. We adapted the method
derived in Gendzwill & Stauffer (1981) to project ellipsoids onto the
celestial sphere.

We denote the 3 ellipsoidal axis lengths as $\lambda_{j}$, $j=1,2,3$. We then
define the diagonal matrix ${\bf\Gamma}$ such that
$\Gamma_{ij}=\delta_{ij}\lambda_{j}^{-2}$, where $\delta$ is a Kronecker
delta. We normalize the corresponding axis directions $\vec{s}_{j}$ and
organize them as rows of a matrix ${\bf S}$, so that $S_{ij}$ is the $j^{\rm
th}$ component of the $i^{\rm th}$ vector. We are projecting along the
$\hat{x}$ unit vector direction, here denoted as component 1, onto the
$\hat{y}-\hat{z}$ plane.

We define the column vector $\vec{m}$ as

$\displaystyle\vec{m}=\left(\hat{x}^{T}{\bf S}^{T}{\bf\Gamma
S}\hat{x}\right)^{-1}\hat{x}^{T}{\bf S}^{T}{\bf\Gamma S}$ (27)

where the pre-factor adopts the normalization that $\vec{m}\cdot\hat{x}=1$. We
then compute vectors $\vec{u}$ and $\vec{v}$ with elements
$u_{j}=\hat{y}\cdot(\vec{m}\times\vec{s}_{j})$ and
$v_{j}=\hat{z}\cdot(\vec{m}\times\vec{s}_{j})$, written alternatively as

$\displaystyle u_{j}$ $\displaystyle=$ $\displaystyle m_{1}S_{j3}-m_{3}S_{j1}$
(28) $\displaystyle v_{j}$ $\displaystyle=$ $\displaystyle
m_{1}S_{j2}-m_{2}S_{j1}.$ (29)

We use these to compute the scalars $A=\vec{u}^{T}{\bf\Gamma}\vec{u}$,
$B=\vec{u}^{T}{\bf\Gamma}\vec{v}$, and $C=\vec{v}^{T}{\bf\Gamma}\vec{v}$.
The orientation angle of the projected ellipse’s primary axis, measured in the
$+\hat{y}$ direction from $\hat{z}$ is

$\displaystyle\tan 2\theta=\frac{-2B}{A-C}$ (30)

And the minor and major axis lengths of the ellipse, $b$ and $a$ are given as:

$\displaystyle\frac{1}{a^{2}}=\frac{A+C}{2}+\frac{A-C}{2\cos{2\theta}}$ (31)
$\displaystyle\frac{1}{b^{2}}=A+C-\frac{1}{a^{2}}$ (32)

To project the shapes on the sky, we rotated the original ellipsoid
eigenvectors using the object’s right ascension and declination, so that
$\hat{x}$ lay along the LOS. This results in the axis lengths and orientation
angle measured East of North for each halo. A function which performs these
operations is available here: github.com/cmlamman/ellipse_alignment

##  Appendix B Expanded Derivations

###  B.1 Shape-Density Correlation

Starting from Equation 16, we can continue the computation as:

$\displaystyle w_{\times}(R)$ $\displaystyle=$
$\displaystyle\frac{\tau}{L}\int dz\
\int\frac{d^{3}q}{(2\pi)^{3}}\left[\frac{{q^{2}/3}-q^{2}_{y}}{q^{2}}\right]\left.e^{-i\vec{q}\cdot\vec{x}}\right|_{\vec{x}=0}\int\frac{d^{3}k}{(2\pi)^{3}}\left.e^{i\vec{k}\cdot\vec{r}}\right|_{\vec{r}=(R,0,z)}\left<\tilde{\rho}^{*}(\vec{q})\tilde{\rho}(\vec{k})\right>$
(33) $\displaystyle=$ $\displaystyle\frac{\tau}{L}\int dz\
\int\frac{d^{3}k}{(2\pi)^{3}}\left[\frac{k^{2}}{3}-k^{2}_{y}\right]k^{-2}P(k)\left.e^{i\vec{k}\cdot\vec{r}}\right|_{\vec{r}=(R,0,z)}.$
(34)

Next, the integral over $z$ creates $\int dz\
\exp(ik_{z}z)=(2\pi)\delta^{D}(k_{z})$. We denote the space of $(k_{x},k_{y})$
as $\vec{K}$, and similarly $\vec{R}$ as $(x,y)$. So we have

$w_{\times}(R)=\frac{\tau}{L}\int\frac{d^{2}K}{(2\pi)^{2}}\left(\frac{K^{2}}{3}-K_{y}^{2}\right)K^{-2}P(K)e^{iK_{x}R}.$
(35)

To simplify this, we introduce

$\Phi(\vec{R})=\int\frac{d^{2}K}{(2\pi)^{2}}\frac{P(K)}{K^{2}}e^{i\vec{K}\cdot\vec{R}},$
(36)

which in turn implies

$w_{\times}(R)=\frac{\tau}{3L}\left(2\partial^{2}_{y}-\partial^{2}_{x}\right)\left.\Phi(\vec{R})\right|_{\vec{R}=R\hat{x}}.$
(37)

$\Phi(\vec{R})$ is isotropic, and can be simplified to a Hankel transform

$\Phi(R)=\int\frac{K\ dK}{2\pi}\frac{P(K)}{K^{2}}J_{0}(KR)$ (38)

with $J_{0}$ being the Bessel function. For a general function $f(R)$, we have
$\partial^{2}f/\partial x^{2}=\partial^{2}f/\partial R^{2}$ and
$\partial^{2}f/\partial y^{2}=(1/R)\partial f/\partial R$. So we have

$w_{\times}(R)=\frac{\tau}{3L}\left(\frac{2}{R}\partial_{R}-\partial^{2}_{R}\right)\Phi(R)=\frac{\tau}{3L}R^{2}\frac{d}{dR}\left[\frac{1}{R^{2}}\Psi(R)\right]$
(39)

where we introduce

$\Psi(R)=-\frac{d\Phi}{dR}=\int\frac{K\;dK}{2\pi}\frac{P(K)}{K}J_{1}(KR),$
(40)

using $dJ_{0}(x)/dx=J_{1}(x)$.

###  B.2 Shape - $\xi_{2}$ Correlation

Details of the derivation of Equation 22

Using $L_{2}(\mu)=(3/2)\mu^{2}-(1/2)$,

$q_{z}^{2}-\frac{q^{2}}{3}=q^{2}\left(\mu_{q}^{2}-\frac{1}{3}\right)=\frac{2q^{2}}{3}L_{2}(\mu_{q})$
(41)

for a 3-d vector $\vec{q}$, and $L_{\ell}=\sqrt{4\pi/(2\ell+1)}Y_{\ell 0}$. We
note that

$\frac{T_{zz}}{2}=\frac{1}{3}\int\frac{d^{3}k}{(2\pi)^{3}}L_{2}(\mu_{k})\tilde{\rho}(\vec{k})e^{i\vec{k}\cdot\vec{r}}.$
(42)

Finally, we have the expansion of a plane wave into spherical harmonics and
spherical Bessel functions:

$e^{i\vec{q}\cdot\vec{r}}=4\pi\sum_{\ell m}i^{\ell}j_{\ell}(qr)Y^{*}_{\ell
m}(\hat{q})Y_{\ell m}(\hat{r}).$ (43)

We then compute $\left<\epsilon_{zz}Q(r)\right>$ as

$\left<\epsilon_{zz}Q(r)\right>={5\tau}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{1}{3}L_{2}(\hat{q})\int\frac{d^{2}\hat{r}}{4\pi}L_{2}(\hat{r})\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{q}\cdot\vec{r}}\left<\tilde{\rho}^{*}(\vec{q})\tilde{\rho}(\vec{k})\right>$
(44)

Converting to power, doing the $\vec{k}$ integral, and expanding the plane
wave yields

$\left<\epsilon_{zz}Q(r)\right>=\frac{5\tau}{3}\int\frac{q^{2}dq}{2\pi^{2}}P(q)\int\frac{d^{2}\hat{q}}{4\pi}L_{2}(\hat{q})\int\frac{d^{2}\hat{r}}{4\pi}L_{2}(\hat{r})4\pi\sum_{\ell
m}i^{\ell}j_{\ell}(qr)Y^{*}_{\ell m}(\hat{q})Y_{\ell m}(\hat{r}).$ (45)

We then can do the two angular integrals, yielding the simpler form:

$\left<\epsilon_{zz}Q(r)\right>=-\frac{\tau}{3}\int\frac{q^{2}dq}{2\pi^{2}}P(q)j_{2}(qr).$
(46)


# A stochastic least action principle applied in the description of black swan
events

Tatiana Cardoso e Bufalo
Department of Physics
Federal University of Lavras
Lavras, MG, Brazil
<EMAIL_ADDRESS>
&R. Bufalo
Department of Physics
Federal University of Lavras
Lavras, MG, Brazil
<EMAIL_ADDRESS>
&Lucas P. G. de Figueiredo
Department of Physics
Federal University of Lavras
Lavras, MG, Brazil
<EMAIL_ADDRESS>
&Qiuping A. Wang
Laboratoire SCIQ, EISEA Numérique et Société,
ESIEA Group
Paris, France
Laboratoire IMMM, CNRS
Le Mans University
Le Mans, France
<EMAIL_ADDRESS>
&Fabio Lucio Alves
Nanjing University
Nanjing, China
<EMAIL_ADDRESS>

###### Abstract

In this paper we present a formulation of the stochastic least action
principle (SAP) to encompass random movements with black swan events (of non
dissipative systems) in terms of heavy tailed distributions. The black swan
events are rare and drastic events, such as earthquakes and financial crisis.
It has been observed that the Tsallis entropy suits well in the description of
the black swan events rather than the Shannon-Boltzman-Gibbs entropy, which is
intrinsically related to the fact that black swan events of physical systems
are proportional to non-local correlations. As a consequence, we could assess
the validity of the path probability distribution obtained using the non-
additive Tsallis entropy.

_Keywords_ Stochastic processes; Black swan events; Stable distributions;
Tsallis q-entropy.

##  1 Introduction

The least action principle (LAP), also known as Hamilton’s principle, is one
of the most important principles in physics. It is related to a fundamental
rule of nature that a system evolves in such a way to extremize an energy
related quantity called action [1]. It is known that LAP applies only to
Hamiltonian systems having energy conservation in deterministic motion [1].
More recently, LAP has been extended to a random motion of Hamiltonian systems
which statistically conserve energy (see ref. [2] and references therein for a
review), a type of stochastic process that can be used to describe diffusion
[3]. This extension has been called stochastic action principle (SAP), which
is motivated by the century old conflict between the second law of
thermodynamics and the Liouville’s theorem in the classical mechanics, as
discussed in the refs. [4, 5], as well as by the question whether and how the
path probability depends on the action in random motion of non-dissipative
systems [6]. The predictions of SAP has been verified by the outcomes of the
numerical experiments of random motion with a Gaussian noise [2, 6], which
confirm, from a probabilistic point of view, the existence of a classical
analogue of the Feynman factor $e^{iA}$ for the path integral formalism [7,
8]. This factor also unveiled a close relationship between the path integral
formalism and the action functional A in classical mechanics.

This work aims to generalize the previous mentioned finding to a different
class of motion, those not described by a Gaussian distribution, since the SAP
formulation is not related to any specific type of random motion. In
particular, we shall focus our attention on examining the validity of the SAP
in the description of black swan events, a type of rare and extreme event [9].
In order to carry out this analysis, we shall consider random noises generated
by a stable distribution, a distribution which can be heavy tailed depending
on the choice of its free parameter (see equation (14)). Essentially, this
modification, corresponds to an anomalous diffusion processes, presenting long
range interactions [10]. A natural framework to describe this kind of non-
local (long range) correlation is the Tsallis $q$-entropy [11].

The Tsallis $q$-entropy was formulated in the context of multifractals and was
initially presented aiming the possibility of generalizing the Shannon-
Boltzmann-Gibbs entropy [11]. Among the many interesting properties that the
use of the Tsallis entropy presents, one can observe that it describes systems
with long range interactions, while the Shannon-Boltzmann-Gibbs entropy is
limited to systems with short-range interactions [12, 13]. This last statement
enunciates that the Tsallis $q$-entropy is the legitimate framework for the
description of anomalous diffusion processes [14]. Moreover, the Tsallis
$q$-entropy has been applied to a wide and vast class of several physical
systems, ranging from solar winds measurements, chaos, financial markets,
quantum information, high energy physics and even to black hole entropy
analysis [15, 16, 17, 18, 19].

Another important point worth to mention about the Tsallis $q$-entropy is that
for suitable values of the parameter $q$, one expects that small probabilities
become more significant in the characterization of the system. As can be
noticed, that is an essential feature to consider in the description of black
swan events. Hence, the properties discussed above strongly support the use of
the Tsallis entropy for the description of the black swan events within the
SAP formalism.

In this paper we present an analysis of the black swan events (of non
dissipative systems) within the stochastic least action principle,
supplemented by the Tsallis $q$-entropy in order to take into account the non-
local (long range) correlations present in this class of physical system. In
section 2, we review the main aspects of the stochastic least action
principle, as well as the presence of constraints, to obtain a path
probability distribution in terms of the maximization of the entropy. In
section 3, we develop the path probability distribution to the random motion,
governed by the Gaussian noise, using the Shannon entropy and the Fokker-
Planck diffusion equation. Furthermore, we verify the validity of the results
by means of computational simulations for the case of null (free particle),
linear (constant force) and quadratic (simple harmonic oscillator - SHO)
potentials. The section 4 presents the main analysis of this work, where we
develop the SAP for the case of black swan events within the Tsallis
$q$-entropy in terms of stable distributions. The results and perspectives are
summarized in the section 5.

##  2 The stochastic least action principle

The classical least action principle states that the action path must be
minimal: $\delta A=0$. However, in a scholastic context, which is
probabilistic, we have to state analogously that [2]

$\overline{\delta A}=\sum_{k}p_{k}\delta A_{k}=0,$ (1)

where $p_{k}$ is the probability distribution chosen to describe the physical
system. Moreover, using the mean action of the paths, i.e.,
$\overline{A}=\sum\limits_{k}p_{k}A_{k}$, and performing a variation of it, we
find

$\delta\overline{A}=\delta Q+\overline{\delta A}.$ (2)

The quantity $\delta Q$ is identified as the varentropy of the system. This
notion comes from the thermodynamics, in which the quantities of a isolated
system as energy ($E$), temperature ($T$) and entropy ($S$) satisfy the
following equation of state

$T\delta S=\delta E=\sum_{i}\delta
p_{i}E_{i}=\delta\overline{E}-\overline{\delta E}.$ (3)

Hence, by comparing the equations (2) and (3), we can establish that

$\delta Q=\frac{\delta S}{\gamma}.$ (4)

Finally, replacing the equations (1) and (4) into equation (2), it results in

$\delta\left(S-\gamma\sum_{k}p_{k}A_{k}+\alpha\sum_{k}p_{k}\right)=0,$ (5)

where we have added the natural constraint, $\sum\limits_{k}p_{k}=1$, by means
of obtaining a normalized probability distribution. In the expression (5), the
quantities $\gamma$ and $\alpha$ are interpreted as Lagrange multipliers.
Therefore, we can conclude that, in order to obtain a path probability
distribution in terms of the stochastic least action principle, we need to
maximize the entropy $S$ subjected to the constraints in (5). In other words,
this statement is the Jaynes maximum entropy principle [20].

##  3 The path probability distribution of random motion with Gaussian noise

One of the most important physical formulation of entropy, which allowed to
establish a relation between the statistical mechanics and the thermodynamics,
is the Shannon-Boltzmann-Gibbs entropy [21]:

$S=-\sum_{k}p_{k}\ln p_{k}.$ (6)

Its physical applications are very abundant, as the information theory [22],
hydrology [23], image processing [24], quantum information [25] and so on. We
can maximize the entropy equation in (6) using the Lagrange formalism. Thus,
considering the constraints (5), we can write the following Lagrangian

$\mathcal{L}=-\sum_{k}p_{k}\ln p_{k}-\gamma\left(\sum_{k}p_{k}A_{k}-\langle
A\rangle\right)+(1-\alpha)\left(\sum_{k}p_{k}-1\right).$ (7)

We can evaluate the Euler-Lagrange equation of motion as

$-\ln p_{i}-\alpha-\gamma\langle A\rangle=0.$ (8)

The equation (8) can be reorganized using the constraints equations, resulting
in

$p_{k}=\frac{1}{Z}e^{-\gamma A_{k}},$ (9)

where $Z=\sum\limits_{i}e^{-\gamma A_{i}}$ is the normalization factor, which
corresponds to a distribution with the shape of the Boltzmann-Gibbs
distribution.

In order to examine the physical content of the equation (9), we perform
computational simulations using the Langevin description of random motion of
particles111A type of mechanical stochastic processes where we can consider
the trajectories as markovian chains [3] and based on that compute its path
probability., a microscopical approach of diffusion processes [3]. This
formalism is based on the solutions of the following equation of motion

$m\frac{d^{2}x}{dt^{2}}=-\frac{dV(x)}{dx}-m\zeta\frac{dx}{dt}-R,$ (10)

where $m$ is the mass of the particles, $x(t)$ is the position, $t$ is the
time, $V$ is the potential, $\zeta$ is the dissipation constant and $R$ is a
random force. We consider a non-dissipative scenario ($\zeta\rightarrow 0$)
where the system can be described by the conservative forces and the Langevin
equation becomes the Newtonian second law of motion with a noise
superposition. The Langevin approach is equivalent to the macroscopic
description using the diffusion equation

$\frac{\partial p}{\partial t}-D\frac{\partial^{2}p}{\partial x^{2}}=0,$ (11)

also known as the Fokker-Planck equation, where $p(x)$ is the position
probability. The solutions of the equation (10),can be written as a
superposition of the classical path $f$ with a random noise $\chi$, resulting
in

$x_{i}=x_{i-1}+\chi_{i}+f(t_{i})-f(t_{i-1}).$ (12)

The explicit expression of the noise $\chi$ is defined in terms of the type of
the diffusion process we are observing. The most known diffusion process in
the literature is the Brownian motion [26], the normal diffusion process
characterized by the Gaussian distribution, i.e., $\chi$ is a random noise
gaussianally distributed.

The simulation is prepared considering the following parameters: the
trajectories are discrete with $10$ time steps with size $dt=10^{-5}s$, the
particles are identical with mass $m=1.39\times 10^{-15}kg$ and the Gaussian
distribution has standard deviation $\sigma=\sqrt{2Ddt}\approx 3\times
10^{-9}m$, where $D=4.3\times 10^{-13}m^{2}/s$ is the diffusion constant. The
chosen values of these parameters correspond to the scenario of particles of
the molecule $SiO_{2}$ in water at normal temperature and pressure conditions,
as suggested in ref. [2]. Under these considerations, we create $10^{3}$
trajectories with width $\delta=\sigma$, launching $10^{6}$ particles for
$10^{2}$ times. The action of the paths is written in terms of a discrete
approximation as

$A_{k}=\sum_{i=1}^{10}\left[\frac{1}{2}m\left(\frac{x_{i+1}-x_{i}}{dt}\right)^{2}-V\left(\frac{x_{i}+x_{i+1}}{2}\right)\right]dt.$
(13)

We can understand the probability of each path as the ratio between the number
of particles that moved along each trajectory and the number of particles that
arrive at the final point. As we can observe in figure 1, the paths’
probabilities behave as predicted by the expression (9) for different
potentials: null ($V=0$), linear ($V=mgx$, with $g=10m/s^{2}$) and quadratic
($V=m\omega^{2}x^{2}/2$, with $\omega=2,5\times 10^{4}rad/s$ [2]).

(a)

(b)

(c)

Figure 1: Simulations of the path probability, logarithm of the probability
versus the action of the paths, with Gaussian noise for three different
scenarios: (a) null potential, (b) linear potential and (c) quadratic
potential.

Furthermore, for the systems analyzed in the simulations, as depicted in
figure 1, we observe that the $\ln P(A)$ and the $A$ are strongly correlated,
verifying therefore the exponential probability distribution in equation (9).
The values of the correlation and the parameter $\gamma$ (Lagrange multipler)
are presented in the table 1. In addition, it was found that by using an
equivalent approach to calculate the path probability to the normal diffusion
process, known as the Wiener path integral [27], that
$\gamma=\frac{1}{2mD}=-8.36\times 10^{26}J^{-1}s^{-1}$ [4], which is in good
agreement with the simulation results, see table 1. If we think in a analogous
way as in the thermodynamics where we have $\gamma\propto T^{-1}$ [21], the
physical meaning of $\gamma$ in our analysis is the same. In the
thermodynamics, we have that the temperature increases the energy of the
states that are distributed by the Boltzmann-Gibbs statistics. On the other
hand, in our ensemble, we observe that the diffusion constant can increase the
action of the (allowed) paths of the system.

Potential | Correlation |  $\gamma$ ($J^{-1}s^{-1}$)
---|---|---
Null | 0.987 | $-6.77\times 10^{26}$
Linear | 0.948 | $-6.15\times 10^{26}$
Quadratic | 0.948 | $-6.97\times 10^{26}$
Table 1: Numerical results of the Gaussian simulation.

##  4 The path probability with black swans events

Throughout this section, we are interested to establish the framework of the
stochastic least action principle for the case of black swan events. In order
to reach this goal, we will examine the path probability generated for this
type of event. The proposed analysis consists in modifying the statistical
distribution of the random noise $\chi$ in the equation (12) to a heavy tailed
distribution instead of a Gaussian distribution [28]. This is based on the
fact that for heavy tailed distributions, the probability of extreme events
are higher as heavier are the tails. That being said, we choose to work with
the symmetric stable distribution due to its close relation with the Gaussian
distribution [10], as can be seen in figure 2 for the sake of illustration.
The characteristic function of the stable distribution is

$\varphi(k)=e^{-(ck)^{\alpha}},$ (14)

where the parameter $\alpha$ is the stability parameter, which controls the
heaviness of the tails and varies within the interval $(0,2]$, such that the
smaller the parameter $\alpha$ is, the heavier are the tails. Moreover, we are
able to recover the Gaussian distribution by considering the following
choices: $\alpha=2$ and $c=\sqrt{\sigma^{2}/2}$, where $\sigma$ corresponds to
the standard deviation. Therefore, we will use the relation between $c$ and
$\sigma$, $c=\sqrt{\sigma^{2}/2}$, to obtain the value of the $c$ parameter of
the stable distribution (14) used on the simulations of black swan events.

Figure 2: Stable distributions with different values of $\alpha$

Although the use of a heavy tailed symmetric stable distribution to describe
black swan events might seems to be naive and straightforward, we observe
that, when performing the simulations with this new random noise distribution,
the path probability distribution does not have a linear correlation with the
action on a log-normal plot. This can be understood because, besides the
change in the random noise probability distribution, we necessarily modify the
diffusion behavior of the physical system. Actually, macroscopically, this
modification means that the diffusion process is anomalous and no longer
described by the ordinary diffusion equation (11). Rather, this diffusion
process is governed by a fractional Fokker-Planck equation,

$\frac{\partial p}{\partial t}-D_{\alpha}\frac{\partial^{\alpha}p}{\partial
x^{\alpha}}=0,$ (15)

where $p(x)$ is the position probability distribution (which is the stable
distribution), $\alpha$ is the parameter defined in the equation (14) and
$D_{\alpha}$ is the related diffusion constant. Since $\alpha$ assumes values
within the interval $(0,2]$, we observe the presence of a fractional
differential operator in the spatial derivative in (15), which is in fact a
non-local operator [29]. This implies the presence of non-local interactions
in our system, which means that the position probability $p$ is affected not
only by its neighborhood, but also by distant points. This kind of phenomena
is observed, for instance, in turbulent diffusion processes [10, 30].

Due to the presence of the non-locality, the path distribution analysis is no
longer valid since the description developed in sections 2 and 3 works well to
short-range interactions [12]. This is a well-known condition in the context
of the Boltzmann-Gibbs distribution to the thermodynamics, but also a
limitation of this formalism. Owing to this drawback, new approaches for
defining entropy have been developed to encompass systems with long-range
correlations [13]. Within this context, a suitable framework is the Tsallis
entropy, also known as $q$-entropy [11, 15]. The main aspect of this proposal
is the presence of an $q$ exponent (which characterizes the physical system in
consideration) in the probability distribution of the system [11] (the
Boltzmann-Gibbs distribution is recovered in the limit $q\to 1$).

Another feature to emphasize about the Tsallis $q$-entropy, extremely
important in the description of the black swan events, is that for suitable
values of the parameter $q$, small probabilities become more significant and
relevant in the description of the system. Furthermore, one can mention that
the Tsallis entropy plays a role in the context of the anomalous diffusion
processes in an equivalent way as the Shannon entropy to the normal diffusion
[14]. Therefore, based on the aforementioned properties, one can say that the
Tsallis entropy is the natural framework to describe black swan events within
the SAP formalism in terms of stable distribution.

The Tsallis entropy can be defined in terms of [15]

$S_{q}=-\sum_{i}p_{i}^{q}\ln_{q}p_{i},$ (16)

where $\ln_{q}x$ is a generalization of the logarithm function

$\ln_{q}x=\frac{x^{1-q}-1}{1-q}.$ (17)

The original logarithm function, and consequently the Shannon entropy, is
recovered in the limit $q\rightarrow 1$. In addition, the inverse of this
function is a generalization of the exponential function

$e_{q}^{x}=\left[1+(1-q)x\right]^{\frac{1}{1-q}}.$ (18)

All the properties of the Shannon-Boltzmann-Gibbs entropy are satisfied by the
Tsallis entropy (16), except the property of additivity. Actually, for two
independent systems $A$ and $B$, the Tsallis entropy has a pseudo additivity
property

$S_{q}(A+B)=S_{q}(A)+S_{q}(B)+(1-q)S_{q}(A)S_{q}(B).$ (19)

Hence, in order to obtain the path probability distribution for the random
motion encompassing black swan events, within the framework developed in
sections 2 and 3, we shall maximize the Tsallis entropy $S_{q}$ equation (16)
subject to the constraints obtained by the stochastic least action principle
(5). Nonetheless, when the Tsallis entropy is considered, the constraints (5)
must be modified to avoid divergences in the partition function (as well as in
the energy of the system) [31]. Moreover, this modification is also important
to the maintenance of the thermodynamic Legendre structure. Hence, under these
considerations, we find the new constraint equation

$\delta\left(S_{q}-\gamma^{\prime}\frac{\sum\limits_{k}p_{k}^{q}A_{k}}{\sum\limits_{k}p_{k}^{q}}+\alpha\sum_{k}p_{k}\right)=0.$
(20)

Furthermore, considering the equation (20) and the Lagrange formalism, the
path probability distribution can be extracted from the Lagrangian

$\mathcal{L}=\left(\frac{1-\sum\limits_{k}p_{k}^{q}}{q-1}\right)+\alpha\left(\sum_{k}p_{k}-1\right)-\alpha\gamma^{\prime}(q-1)\left(\sum_{k}p_{k}A_{k}-\langle
A\rangle\right).$ (21)

Thus, evaluating the Euler-Lagrange equations of motion, we determine the
following distribution

$p_{k}=\frac{e_{q}^{-\gamma A_{k}}}{\sum\limits_{i}e_{q}^{-\gamma A_{i}}}$
(22)

where $\gamma\equiv\frac{\gamma^{\prime}}{\sum\limits_{k}p_{k}^{q}}$. We have
also performed the replacement $A_{k}\to A_{k}-\langle A\rangle$, being
$\langle A\rangle=A_{\mbox{clas}}$, since the noise distribution is symmetric.

Therefore, once we have rewritten the path probability distribution, now
considering the Tsallis entropy, we can perform the computational simulations
considering the potentials: null, linear and quadratic, for different values
of $\alpha$. This analysis will allow us to verify the validity of the path
probability distribution (22) as shown below. In addition, since we are
interested on the rare events that lie within the heavy tailed distribution,
we have increased the number of paths from $10^{3}$ to $10^{5}$. The results
from the computational simulations for the logarithm of the probability ($\ln
P$) versus the action of the paths ($A$) can be seen in the figures 3, 4 and
5.

(a) $\alpha=2$

(b) $\alpha=1.75$

(c) $\alpha=1.5$

(d) $\alpha=1.25$

(e) $\alpha=1$

(f) $\alpha=0.75$

(g) $\alpha=0.5$

(h) $\alpha=0.25$

Figure 3: Computational simulations for the free particle scenario considering
the stable distribution for different values of $\alpha$.

(a) $\alpha=2$

(b) $\alpha=1.75$

(c) $\alpha=1.5$

(d) $\alpha=1.25$

(e) $\alpha=1$

(f) $\alpha=0.75$

(g) $\alpha=0.5$

(h) $\alpha=0.25$

Figure 4: Computational simulations for a particle under a linear potential
scenario considering the stable distribution for different values of $\alpha$.

(a) $\alpha=2$

(b) $\alpha=1.75$

(c) $\alpha=1.5$

(d) $\alpha=1.25$

(e) $\alpha=1$

(f) $\alpha=0.75$

(g) $\alpha=0.5$

(h) $\alpha=0.25$

Figure 5: Computational simulations for a particle under a quadratic potential scenario considering the stable distribution for different values of $\alpha$. | Null | Linear | Quadratic
---|---|---|---
$\alpha$ | q | q-error | q | q-error | q | q-error
2 | 0.984 | 0.002 | 0.963 | 0.003 | 1.009 | 0.004
1.75 | 1.113 | 0.002 | 1.092 | 0.003 | 1.124 | 0.004
1.5 | 1.227 | 0.003 | 1.203 | 0.004 | 1.236 | 0.004
1.25 | 1.315 | 0.005 | 1.286 | 0.005 | 1.322 | 0.006
1 | 1.396 | 0.007 | 1.358 | 0.007 | 1.398 | 0.008
0.75 | 1.497 | 0.012 | 1.448 | 0.013 | 1.492 | 0.013
0.5 | 1.536 | 0.027 | 1.558 | 0.035 | 1.593 | 0.032
0.25 | 1.749 | 0.078 | 1.633 | 0.117 | 1.695 | 0.072
Table 2: Values of $q$ and its uncertainty for the null, linear and quadratic
potential scenarios and different values of the parameter $\alpha$.

The numerical results of the parameter $q$ are summarized in the table 2 for
different values of $\alpha$ and the cases of null, linear and quadratic
potentials. The outcomes of our analysis allow us to strongly conclude that
the path probability of particles under random motion with stable
distributions are, in fact, distributed by the maximization of the Tsallis
entropy subjected to the constraints of the stochastic least action principle.
In other words, it means that in order to describe systems with drastic and
rare phenomena (named as black swan events), using stable distributions and
characterized by the action, the Tsallis entropy is preferred instead of the
Shannon entropy, since we are dealing with non-local correlations.

Furthermore, from the results presented in the table 2, it is possible to
observe a relation between the parameters $\alpha$ and $q$, which behaves with
good agreement with the data as $q=\sqrt{3-\alpha}$, see figure 6. Actually,
this result can be presented as the following statement: the heavier the tails
are, stronger is the non-locality effect. This conclusion is rather expected,
since $1-q$ can be understood as a measure of how much the system “diverges”
from the (normal) behavior, once in the limit $q\rightarrow 1$, the Shannon
entropy is recovered.

Figure 6: The dots represent the values of $q$ and $\alpha$ to the simulations
with different potentials: blue dot for null potential, orange dot for linear
potential and green dot for quadratic potential. The continuous line
represents the approximated behavior of these two parameters as
$q=\sqrt{3-\alpha}$.

##  5 Conclusions

In this work we have generalized the stochastic least action principle to
describe the so-called black swan events (of non dissipative systems) with the
use of the Tsallis $q$-entropy and heavy tailed distribution. This approach
led to a new path probability distribution, in terms of the $q$-exponential,
for the random motion of particles. Computational simulations have shown the
agreement of the present formalism with the path probability of random motion
with stable distribution for the scenarios of null, linear and quadratic
potential. These results support the choice of using the Tsallis entropy,
since we are dealing with long range interactions of the anomalous diffusion.
Additionally, the relation between the stable distribution and the Tsallis
entropy seems to be well modeled with the form $q=\sqrt{3-\alpha}$, which can
be interpreted as an equivalence between the long range interactions and the
occurrence of the black swan events. Therefore, we can state that within the
SAP formalism, the Tsallis $q$-entropy provides the natural environment for
the description of the black swan events.

It is worth recalling that the SAP’s main strength consists in prescribing the
path probability for different types of energy conserving random motion. Its
application to thermodynamic systems can help to shed light on the disparities
between the thermodynamics laws and classical mechanics laws [4, 5]. Here, we
have only investigated black swan events in energy-conserving dynamical
systems. Naturally, a direct extension would be a formulation of black swan
events within dissipative systems with the purpose of treating random systems
in a wider framework [6, 32]. We can mention some scenarios where SAP can be
suitably applied to with the technique of path integrals [8]: financial Black-
Scholes formalism [33], models of biological evolution [34] and the human body
neural system [35], among others. Efforts in these directions are now under
development and will be reported elsewhere.

### Acknowledgments

This study was funded partially by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001. R.B.
acknowledges partial support from Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq Projects No. 305427/2019-9 and No.
421886/2018-8) and Fundação de Amparo à Pesquisa do Estado de Minas Gerais
(FAPEMIG Project No. APQ-01142-17).

## References

  * [1] Herbert Goldstein, Charles Poole, and John Safko.  Classical mechanics. American Association of Physics Teachers, 2002.
  * [2] Tongling Lin.  Path probability and an extension of least action principle to random motion. PhD thesis, Le Mans, 2013.
  * [3] Masud Chaichian and Andrei Demichev.  Path integrals in physics: Volume I stochastic processes and quantum mechanics. CRC Press, 2018.
  * [4] Qiuping A. Wang.  Maximum path information and the principle of least action for chaotic system. Chaos, Solitons & Fractals, 23(4):1253–1258, 2005.
  * [5] Qiuping A. Wang and Aziz El Kaabouchiu.  “From random motion of hamiltonian systems to boltzmann’s h theorem and second law of thermodynamics: a pathway by path probability.”  Entropy, 16(2):885–894, 2014.
  * [6] T.L. Lin, R Wang, W.P. Bi, A. El Kaabouchi, C. Pujos, F. Calvayrac, and Q. Alexandre Wang.  “ Path probability distribution of stochastic motion of non dissipative systems: a classical analog of feynman factor of path integral.”  Chaos, Solitons & Fractals, 57:129–136, 2013.
  * [7] Richard P. Feynman, Albert R. Hibbs, and Daniel F. Styer.  Quantum mechanics and path integrals. Courier Corporation, 2010.
  * [8] Hagen Kleinert.  Path integrals in quantum mechanics, statistics, polymer physics, and financial markets. World scientific, 2009.
  * [9] Nassim Nicholas Taleb.  The black swan: The impact of the highly improbable, volume 2. Random house, 2007.
  * [10] R. Klages, G. Radons, and I.M. Sokolov.  Anomalous Transport: Foundations and Applications. Wiley, 2008.
  * [11] Constantino Tsallis.  “Possible generalization of boltzmann-gibbs statistics.”  Journal of statistical physics, 52(1):479–487, 1988.
  * [12] Peter T. Landsberg.  Thermodynamics and statistical mechanics. Oxford University Press, 1978.
  * [13] Sumiyoshi Abe and Yuko Okamoto.  Nonextensive statistical mechanics and its applications. Springer Science & Business Media, 2001.
  * [14] C. Tsallis and E.K. Lenzi.  “Anomalous diffusion: nonlinear fractional Fokker–Planck equation.”  Chemical Physics, 284(1-2):341–347, 2002.
  * [15] Murray Gell-Mann and Constantino Tsallis.  Nonextensive entropy: interdisciplinary applications. Oxford University Press, 2004.
  * [16] Constantino Tsallis.  “Mecânica estatística de sistemas complexos.”  Revista Brasileira de Ensino de Física, 43, 2021.
  * [17] Anastasios Bountis, J.J.P. Veerman, and Franco Vivaldi.  “Cauchy distributions for the integrable standard map.”  Physics Letters A, 384(26):126659, 2020.
  * [18] Donghua Xiao, Xingkun Peng, Youlong Yuan, Qilong Cai, Huibin Qiu, Tianyi Hu, Honghui Zhang, Shengfa Wu, Xiaobin Li, Jinming Chang, et al.  “Innovation for measuring the distribution function with nonextensive single electric probe.”  AIP Advances, 11(8):085228, 2021.
  * [19] Abhishek Majhi.  “ Non-extensive statistical mechanics and black hole entropy from quantum geometry. ”  Physics Letters B, 775:32–36, 2017.
  * [20] Edwin T. Jaynes.  “Information theory and statistical mechanics.”  Physical review, 106(4):620, 1957.
  * [21] Frederick Reif.  Fundamentals of statistical and thermal physics. Waveland Press, 2009.
  * [22] Claude Elwood Shannon.  “ A mathematical theory of communication.”  The Bell system technical journal, 27(3):379–423, 1948.
  * [23] Jagat Narain Kapur and Hiremaglur K. Kesavan.  “Entropy optimization principles and their applications.”  In Entropy and energy dissipation in water resources, pages 3–20. Springer, 1992.
  * [24] Stephen F. Gull and John Skilling.  “Maximum entropy method in image processing.”  In IEE Proceedings F-Communications, Radar and Signal Processing, volume 131, pages 646–659. IET, 1984.
  * [25] Vlatko Vedral.  “ The role of relative entropy in quantum information theory.”  Reviews of Modern Physics, 74(1):197, 2002.
  * [26] A. Einstein.  “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen.”  Annalen der Physik, 322:549–560, 1905.
  * [27] Robert M. Mazo.  Brownian motion: fluctuations, dynamics, and applications. Number 112. Oxford University Press on Demand, 2002.
  * [28] Hagen Kleinert.  “Quantum field theory of black-swan events.”  Foundations of physics, 44(5):546–556, 2014.
  * [29] Rudolf Hilfer.  Applications of fractional calculus in physics. World scientific, 2000.
  * [30] Paolo Paradisi, Rita Cesari, Francesco Mainardi, and Francesco Tampieri.  “The fractional fick’s law for non-local transport processes.”  Physica A: Statistical Mechanics and its Applications, 293(1-2):130–142, 2001.
  * [31] Constantino Tsallis, Renio S. Mendes, and Anel R. Plastino.  “The role of constraints within generalized nonextensive statistics.”  Physica A: Statistical Mechanics and its Applications, 261(3-4):534–554, 1998.
  * [32] Aziz El Kaabouchi and Qiuping A. Wang.  “Least action principle and stochastic motion: a generic derivation of path probability.”  In Journal of Physics. Conference Series (Online), volume 604, 2015\.
  * [33] Fischer Black and Myron Scholes.  “The pricing of options and corporate liabilities.”  Journal of political economy, 81(3):637–654, 1973.
  * [34] Alexander Pease.  QIB: Connecting Feynman Path Integrals to Biophysical Problems. 2011\.
  * [35] Vladimir G. Ivancevic and Tijana T. Ivancevic.  Quantum leap: From Dirac and Feynman, across the universe, to human body and mind. World Scientific, 2008.