# Out of sight, out of mind? The impact of correlated clustering in
substructure lensing

Alexandres Lazar1, James S. Bullock1, Michael Boylan-Kolchin2, Robert
Feldmann3, Onur Çatmabacak3, and Leonidas Moustakas4
1Department of Physics and Astronomy, University of California, Irvine, CA
92697 USA
2Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Stop C1400, Austin, Texas 78712-1205, USA
3Institute for Computational Science, University of Zurich, Zurich CH-8057,
Switzerland
4Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Dr., Pasadena CA 91109, USA<EMAIL_ADDRESS>
(Working Draft)

###### Abstract

A promising route for revealing the existence of dark matter structures on
mass scales smaller than the faintest galaxies is through their effect on
strong gravitational lenses. We examine the role of local, lens-proximate
clustering in boosting the lensing probability relative to contributions from
substructure and unclustered line-of-sight (LOS) halos. Using two cosmological
simulations that can resolve halo masses of $M_{\rm halo}\simeq 10^{9}\
M_{\odot}$ (in a simulation box of length $L_{\rm box}{\sim}100\,{\rm Mpc}$)
and $10^{7}\ M_{\odot}$ ($L_{\rm box}\sim 20\,{\rm Mpc}$), we demonstrate that
clustering in the vicinity of the lens host produces a clear enhancement
relative to an assumption of unclustered halos that persists to $>20\,R_{\rm
vir}$. This enhancement exceeds estimates that use a two-halo term to account
for clustering, particularly within $2-5\,R_{\rm vir}$. We provide an analytic
expression for this excess, clustered contribution. We find that local
clustering boosts the expected count of $10^{9}\ M_{\odot}$ perturbing halos
by ${\sim}35\%$ compared to substructure alone, a result that will
significantly enhance expected signals for low-redshift ($z_{l}\simeq 0.2$)
lenses, where substructure contributes substantially compared to LOS halos. We
also find that the orientation of the lens with respect to the line of sight
(e.g., whether the line of sight passes through the major axis of the lens)
can also have a significant effect on the lensing signal, boosting counts by
an additional $\sim 50\%$ compared to a random orientations. This could be
important if discovered lenses are biased to be oriented along their principal
axis.

###### keywords:

cosmology:theory – dark matter – large-scale structure – gravitational
lensing: strong

††pubyear: 2020††pagerange: Out of sight, out of mind? The impact of
correlated clustering in substructure lensing–C

##  1 Introduction

The large-scale clustering of galaxies provides important constraints on the
makeup and evolution of the Universe (e.g. Geller & Huchra, 1989; Bond et al.,
1996; Tegmark et al., 2004; Sánchez et al., 2006). The dark energy plus cold
dark matter (CDM) paradigm, $\Lambda$CDM, is entrenched as the benchmark model
for the theory of galaxy formation based largely on its success in matching
observed large-scale structure. For years, cosmological $N$-body simulations
that incorporate only gravitational dynamics (dark matter only, DMO,
simulations) have served as crucial tools for understanding the $\Lambda$CDM
model, and have been used to understand the detailed clustering of galaxies.
When introducing full galaxy formation physics, cosmological simulations are
able to match observed clustering statistics as a function of galaxy type as
well (e.g. Vogelsberger et al., 2014; Genel et al., 2014; Schaye et al., 2015;
Khandai et al., 2015; Dubois et al., 2016; Dolag et al., 2016; Springel et
al., 2018; Vogelsberger et al., 2020) however some discrepancies on smaller
scales exist and motivate the exploration of alternative models (Bullock &
Boylan-Kolchin, 2017; Meneghetti et al., 2020).

A feature of CDM that profoundly separates it from many other dark matter
models is that CDM predicts a rich abundance of low-mass dark matter halos
$M_{\rm halo}<10^{6}\ M_{\odot}$ (Press & Schechter, 1974; Bullock & Boylan-
Kolchin, 2017). In cosmologies that include warm dark matter (WDM), for
example, the power spectrum is suppressed on scales smaller than a value set
by the WDM particle mass (Bode et al., 2001; Bozek et al., 2016). For a
thermal WDM particle of mass $m_{\rm thm}\lesssim 5\ \rm kev$, the formation
of halos $<10^{7}\ M_{\odot}$ (e.g. Schneider et al., 2013; Horiuchi et al.,
2016) is suppressed. Therefore, if halos below $\simeq 10^{7}\ M_{\odot}$ are
detected, this would impose significant constraints on both the dark matter
power spectrum and the particle nature of dark matter (see Bertone & Tait,
2018).

Dark matter halos of sufficiently low mass are expected to be unable to form
stars or retain baryons in the presence of a cosmological photoionizing
background (e.g. Efstathiou, 1992; Bullock et al., 2000). The detection of
these starless halos, with the abundance and density structure predicted by
simulations, would be triumphant for the CDM model. One way of inferring the
presence of these low-mass objects is by their influence on cold, low-velocity
stellar streams in the Milky Way (Ibata et al., 2002; Carlberg, 2009; Yoon et
al., 2011). Recently Banik et al. (2019) argued that the observed
perturbations of the MW’s stellar streams can only be explained by a
population of subhalos in CDM. They set constraints to alternative dark matter
models for halos down the mass function, notably setting a lower limit on the
mass of warm dark matter thermal relics $m_{\rm thm}\gtrsim[4.6-6.3]\ \rm
keV$. In order to provide tighter constraints for substructure down to $M_{\rm
halo}\simeq 10^{5-6}M_{\odot}$ populating the MW, a larger sample cold streams
would be needed. Another proposed approach to detecting these low-mass halos
could be through the kinematics of stars in the Milky Way’s disk (Feldmann &
Spolyar, 2015).

Currently, the field of strong gravitational lensing offers to be another tool
for the indirect detection of low-mass, starless halos of masses $\simeq
10^{6-8}\ M_{\odot}$ (Dalal & Kochanek, 2002; Koopmans, 2005; Vegetti et al.,
2010, 2014; Li et al., 2016; Nierenberg et al., 2017). Lensing perturbations
can arise from both subhalos within the lens host and from small halos found
outside of the virial radius that perturb the light from source to the
observer (dubbed “line-of-sight” (LOS) halos; Li et al., 2017; Despali et al.,
2018). Notably, the field of substructure lensing offers tantalizing
prospects, as in the near future, we expect both a gross increase in the
number of lenses as well as a boost in resolution for instrument sensitivity.
Forecasts suggest that the Dark Energy Survey (DES), LSST and EUCLID should
discover hundreds of thousands of galaxy-galaxy lensing systems (Collett,
2015). The Nancy Grace Roman Space Telescope (RST) also has the potential of
providing complementary catalogs of lens images (Weiner et al., 2020).
Additionally, the detection of halos with $M_{\rm halo}\sim 10^{7-8}\
M_{\odot}$ might be possible with JWST via quasar flux ratio anomalies
(MacLeod et al., 2013). As of now, reported detections using ALMA reach the
mass scale of classical Milky Way satellites ($M_{\rm halo}\sim
10^{10}M_{\odot}$), with constraints on subhalos an order of magnitude
smaller. In the future, further observations and improved constraints may
significantly improve these limits (Hezaveh et al., 2016) and offer tighter
constraints on the warm dark matter mass (He et al., 2020), especially when
combined with Milky Way satellite constraints (Nadler et al., 2021).

The expected count of subhalos that exist within the virial radius of the lens
system has been studied rigorously. Studies of subhalo populations and their
effect on lensing have previously relied on DMO simulations (e.g. Bradač et
al., 2002; Xu et al., 2009; Xu et al., 2012; Metcalf & Amara, 2012; Vegetti et
al., 2014). More recently, however, simulations that implement full galaxy
formation physics show that small subhalos are actually depleted with respect
to DMO simulations, owing to interactions with the central galaxy (Brooks &
Zolotov, 2014; Wetzel et al., 2016; Zhu et al., 2016; Graus et al., 2018;
Kelley et al., 2019; Richings et al., 2021). Notably, Garrison-Kimmel et al.
(2017) showed that it is the central galaxy potential itself, not feedback,
that drives most of the factor of $\sim 2$ difference in subhalo counts
between DMO and full physics simulations for Milky Way-mass halos ($M_{\rm
vir}\simeq 10^{12}\ M_{\odot}$). For lens-mass halos of interest ($M_{\rm
vir}\simeq 10^{13}\ M_{\odot}$), Despali & Vegetti (2017) used both DMO and
full physics from the EAGLE (Schaye et al., 2015) and Illustris (Vogelsberger
et al., 2014) simulations to investigate predictions for subhalo lensing and
found that simulations with full galaxy formation physics reduces the average
expected substructure counts. The substructure analysis done in Graus et al.
(2018) using Illustris found similar results.

Given the expected depletion in subhalo counts seen in full-physics
simulations, the contribution of lensing signals from LOS halos has been
recognized as ever more important. If LOS halos dominate the signal from a
given lens, then uncertainties are reduced substantially because the
contribution of the LOS halos can be accurately calculated, independent of the
effect of baryonic physics, for a variety of cosmologies. Efforts have been
made to understand the contribution of the LOS structure on the flux-ratio
anomalies in lensed quasars (e.g. Metcalf, 2005; Metcalf & Amara, 2012; Xu et
al., 2012; Inoue & Takahashi, 2012; Inoue et al., 2015; Xu et al., 2015;
Inoue, 2016; Gilman et al., 2018, 2019). It has become increasingly apparent
that LOS halos should dominate the signal for more distant lenses ($z_{l}\sim
0.5$) while the contribution from subhalos should be non-negligible for more
local lenses ($z_{l}\sim 0.2$; Despali & Vegetti 2017).

Our objective in this article is to understand and quantify an effect not
discussed in most previous work: how does local, correlated structure, in the
vicinity of the lens host halo, impact the expected lensing signal? While such
an effect has been discussed before in the context of weak lensing (D’Aloisio
et al., 2014), its impact in strong lensing remains elusive. This effect will
be most important for low-redshift lenses, where subhalos are known to
contribute non-trivially compared to the LOS count. We use a suite
cosmological simulations, including those that include both DMO and and full
galaxy formation physics, to explore this question. Specifically, we quantify
correlated structure in lens-centered projections of targets lens halos of
$M_{\rm vir}\simeq 10^{13}\,M_{\odot}$ at redshift $z=0.2$, corresponding with
the benchmark sample discussed in Vegetti et al. (2014). This is done using
two simulation projects: The first, from the public IllustrisTNG project
(Nelson et al., 2018), includes both DMO and full-physics versions and
resolves halos down to $M_{\rm halo}=10^{9}\,M_{\odot}$ in a fairly large
cosmological volume. The second is a DMO simulation evolved in much smaller
cosmological volume that is able resolve dark matter halos down to $M_{\rm
halo}=10^{7}\,M_{\odot}$.

This paper is structured as follows. Section 2 introduces our set of
simulations, provides a description of the selected sample of halos, and
outlines our methodology for counting structures along projected line-of-
sights in the simulations. Section 3 provides results on structure and
explores how the lens-host orientation can boost the amount of structure
expected along lens-centered projections. Finally, we summarize our results in
Section 4. Note that in strong lensing literature, terms such as “field halo”
and “line-of-sight (LOS) halo” typically refer to the same things. We will use
these terms interchangeably throughout this work to refer to halos having a
volume density equal to the cosmological mean density of halos at that mass.

##  2 Numerical Methodology

Our $\Lambda$CDM predictions rely on two sets of simulations. The primary set
comes from the public catalogs of the IllustrisTNG project111The Illustris
data is publicly available at https://www.tng-project.org/ (Nelson et al.,
2018; Marinacci et al., 2018; Springel et al., 2018; Naiman et al., 2018;
Pillepich et al., 2018). As described in Section 2.1, these simulations allow
us to identify halo populations robustly to masses greater than $10^{9}\
M_{\odot}$ for a large-scale environment both with gravitational physics alone
and with full galaxy formation physics. The second simulation, introduced in
Section 2.2, is a DMO version of a new simulation suite called FIREbox, which
is part of the Feedback in Realistic Environments (FIRE) project.222The FIRE
project website: http://fire.northwestern.edu The mass functions of these
simulations are presented in Fig. 1. Both simulations assume a Planck 2015
cosmology (Planck Collaboration et al., 2016): $\Omega_{\rm m}=0.3089$,
$\Omega_{\Lambda}=0.6911$, $\Omega_{\rm b}=0.0486$, $\sigma_{8}=0.8159$, and
$h=0.6774$.

Our analysis relies on halo catalogs at redshift $z=0.2$, which coincides with
a typical redshift of the lens-host sample explored in Vegetti et al. (2014).
Dark matter halos are defined to be spherical systems with a virial radius,
$R_{\rm vir}$, inside which the density is equal to the average density of
$\Delta_{\rm vir}(z)\rho_{\rm crit}(z)$, where $\Delta_{\rm vir}(z)$ is the
virial overdensity defined by Bryan & Norman (1998) and $\rho_{\rm crit}(z)$
is the critical density of the Universe. The virial mass, $M_{\rm vir}$, is
then the total mass enclosed in a sphere of radius $R_{\rm vir}$.

In what follows, we discuss two types of halos: First, we have massive target
halos, chosen to mimic lens-galaxy hosts: $M_{\rm vir}\in[0.8-2]\times
10^{13}\ M_{\odot}$. Second, we have low-mass perturber halos, which can
either be subhalos or small halos that exist somewhere outside of the host’s
virial radius and along the LOS from the observer. We quote their masses
(taken from the halo catalogs) using the symbol $M_{\rm halo}$, since for
subhalos $M_{\rm vir}$ is not physically relevant.

###  2.1 The IllustrisTNG simulations

The IllustrisTNG (TNG) suite of cosmological simulations was run using the
moving-mesh code Arepo (Springel, 2010). The runs with full galaxy formation
physics use an updated version of the Illustris model (Weinberger et al.,
2017; Pillepich et al., 2018). We use the high resolution set of publicly-
available simulations,333The highest resolution box is actually TNG50-1, but
was not publicly available by the time this article was submitted. TNG100-1,
which has a comoving box of side length $L_{\rm box}^{\rm
TNG}=75\,h^{-1}\mathrm{cMpc}=110.7\,{\rm cMpc}$. The full physics run, TNG100,
has a dark matter particle mass of $m_{\rm dm}=7.5\times 10^{6}\ M_{\odot}$
and a gas particle mass of $m_{\rm gas}=1.4\times 10^{6}\ M_{\odot}$. The
high-resolution simulation that uses dark matter only (DMO) physics, TNG100DM,
has the same box size, but treats baryonic matter as collisionless particles,
in which gives the simulation a particles mass is $m_{\rm dm}=8.9\times
10^{6}\ M_{\odot}$. When comparing between TNG100 and TNG100DM, we account for
the excess baryonic mass in the DMO simulation by introducing a factor $m_{\rm
dm}\rightarrow(1-f_{b})m_{\rm dm}$, where $f_{b}:=\Omega_{\rm b}/\Omega_{\rm
m}$ is the cosmic baryon fraction.

For both TNG100 and TNG100DM, the dark matter halo catalogs were constructed
using a friends-of-friends (FOF) algorithm (Davis et al., 1985) with a linking
length $0.2$ times the mean inter-particle spacing. Gravitationally bound
substructures are identified through the SUBFIND halo finder (Springel et al.,
2001). These subhalos have a dark matter mass that is gravitationally bound to
itself but not to any other subhalos found in the same FOF host or the host
itself. As shown in Fig. 1, the $z=0.2$ mass functions of the TNG100 (magenta)
and TNG100DM (cyan) SUBFIND subhalos match well with the analytical Sheth &
Tormen (2002) mass function (dotted dashed) at that redshift. Both simulations
become incomplete just below $M_{\rm halo}\simeq 10^{8.6}\ M_{\odot}$. In
order to be conservative, we restrict our analysis to subhalos and field halos
with resolved masses $M_{\rm halo}>10^{9}\ M_{\odot}$. These halos contain at
least 134 dark matter particles. TNG100 and TNG100DM have 166 and 168 lens-
target halos, respectively, where in TNG100, these halos typically host
galaxies with $M_{\star}\simeq 10^{11}\ M_{\odot}$, which is consistent with
the benchmark sample in Vegetti et al. (2014).

The public SUBFIND catalogs in TNG100 include several baryon-dominated
“subhalos," many of which contain no dark matter. Most of these baryon-
dominated objects exist within ${\sim}20$ kpc of host galaxies and appear to
be baryonic fragments numerically identified by SUBFIND rather than galaxies
associated with dark matter subhalos. While baryonic clumps could induce
perturbations detectable in lens images (Gilman et al., 2017; Hsueh et al.,
2017; Hsueh et al., 2018; He et al., 2018), this type of object is not the
focus of our analysis; we are interested in the search for low-mass dark
matter structures. For this reason, we exclude systems with a ratio of total
baryonic mass to dark matter mass that is more that twice the cosmic baryon
fraction, ($M_{\rm bar}/M_{\rm halo}>2\times f_{\rm b}$) in the substructure
analysis that follows.

Figure 1: — The abundance of resolved dark matter halos. Shown are the halo
mass functions of the simulations at $z=0.2$. Results from all three
simulation are in agreement with the analytical prediction of Sheth & Tormen
(2002). The arrows pointing to the $x$-axis designate the minimally resolved
mass considered in each simulation, indicated by the color coding.

###  2.2 The FIREbox DMO simulation

FIREBox is a new effort within the FIRE collaboration to simulate cosmological
volumes of $L_{\rm box}^{\rm FIRE}=15\ h^{-1}\rm cMpc=22.14\ cMpc$ at the
resolution of FIRE zoom-in simulations (Feldmann et al., in prep). Our
analysis uses the results from a DMO version of the FIREbox initial
conditions, ran with $2048^{3}$ particles, which we dub FIREboxDM in this
article. FIREboxDM has a dark matter particle mass of $m_{\rm DM}=5.0\times
10^{4}\ M_{\odot}$ and was run with a Plummer softening length of
$\epsilon_{\rm DM}=40\ \rm pc$. Since we will work with a DMO simulation to
complement the analysis done in the two TNG simulations, we again account for
baryonic mass by multiplying a by conversion factor of $(1-f_{\rm b})$, as we
do for TNG100DM.These simulations are complete, conservatively, to halo masses
down to $10^{7}\ M_{\odot}$ (see Fig. 1).

The halo catalog for FIREboxDM was generated using ROCKSTAR (Behroozi et al.,
2013). This method finds dark matter halos using a hierarchically adaptive
refinement of FOF in a 6-dimensional phase-space and one time dimension. We
set our ROCKSTAR halo finding parameters to be comparable to those use in the
TNG catalogs generated by SUBFIND. Specifically, we use a FOF linking length
of $b=0.20$ and include only halos that have at least $100$ dark matter
particles. We also set the criteria for unbound particle fraction rejection to
$70\%$ instead of the default $50\%$, as explored in the Appendix of Graus et
al. (2018). Doing so minimizes ambiguities associated with using different
halo finders for computing halo masses.444We have also experimented with
higher unbound fractions of $90\%$ and $95\%$ with a fixed value of $b=0.20$.
We chose the unbound fraction of $70\%$ and a FOF linking of $b=0.20$ because
they provide the best match with TNG100DM catalogs for $M_{\rm halo}>10^{9}\
M_{\odot}$. With these choices, FIREbox contains three target halos with mass
$M_{\rm vir}\in[0.8-2]\times 10^{13}\ M_{\odot}$.

Returning to Fig. 1, the mass function of FIREboxDM at $z=0.2$ is shown as the
solid gray curve and agrees well with the analytical mass function down to
halo masses with $10^{7-8}\ M_{\odot}$. With this in mind, the use of
FIREboxDM in our analysis will be restricted to two sets of subhalo masses: a
sample of halo masses down to $10^{7}\ M_{\odot}$ (${\sim}10^{2}$ particles)
and a sample of down to $10^{8}\ M_{\odot}$ (${\sim}10^{3}$ particles). For a
more stringent check, Appendix C compares the subhalo mass function and the
subhalo $V_{\rm max}$ function of the three target-lens halos of FIREboxDM to
the same mass target-lens halos in TNG100DM. We find excellent agreement
between substructure statistics between our halo samples.

Figure 2:  A cartoon depiction of the analysis we perform for each lens-
centered cylinder in a simulation box. The lens-mass halo, centered at the
origin of the box, is illustrated as the large pink sphere, while subhalos are
illustrated as smaller halos. With the cylinder held at a fixed $\mathcal{R}$,
the cylinder is varied by increments of $\Delta{\ell}$ ($\Delta\ell/2$ at each
end of the cylinder) until reaching the edge of the box. Note that the radius
of the cylinder has been exaggerated for clarity. In practice, $\mathcal{R}\ll
R_{\rm vir}$.  Figure 3:  — The importance of correlated structure in sub-
galactic lensing. The upper and middle panels depict the side-view of a
cylinder length of $\ell=10\ \rm Mpc$ and radius $\mathcal{R}=500\ \rm kpc$
centered on a lens-mass host halo with radius $R_{\rm vir}\simeq 500\ \rm kpc$
at $z=0.2$ in the TNG100 simulation. The points show locations of $M_{\rm
halo}>10^{9}\ M_{\odot}$ halos within the cylinder, color coded by relative
distance from the host out to the size of the halo $\ell=2R_{\rm vir}$ (or
equivalently to $r=R_{\rm vir}$; cyan points), within $\ell=[2-4]R_{\rm vir}$
($r=[1-2]R_{\rm vir}$; magenta points), and everything else outside
$\ell=4R_{\rm vir}$ out to $\ell=10$ Mpc (black points). The figure beneath it
shows an identical cylinder that samples a representative region of the
simulation box using the same color scheme; the area of each point is
proportional to that halo’s mass. The square plots are the same cylinders
shown face on, centered on the lens (left) and centered randomly (right). The
cylinder radius $\mathcal{R}=500\,{\rm kpc}$ used for this figure is, for
illustrative purposes, much larger than the typical Einstein radius of the
host ($\mathcal{R}\approx 10\,{\rm kpc}$; we use the latter value for our
actual analysis. Local, correlated perturbers around the lens are highly
significant (compare the magenta points in the left verses right panels on the
bottom).  Figure 4: — Structure excess along lens-centered projections: The
mean differential count $\langle dN/d\ell\rangle$ of small halos within
cylinders of increasing length $\ell+\Delta{\ell}$ centered on the lens-mass
host halos ($M_{\rm vir}\simeq 10^{13}M_{\odot}$, solid colored), normalized
by the expected average background counts via Sheth & Tormen (2002) halo mass
function (dotted black). The cylinders mirror those shown in Fig 3, but now
with a radius comparable to the Einstein radius of the lens-mass halos,
$\mathcal{R}=10\ \rm kpc$ (and varying length $\ell$). The gray vertical
dashed line marks the typical outer region of the halo $(\ell=2R_{\rm
vir}\leftrightarrow r=R_{\rm vir})$. Inside the lens-mass host halos
($\ell<2R_{\rm vir}$), the differential counts are self-similar for the DMO
simulations. For TNG100 (solid magenta), the counts decrease by almost a
factor of 10 because of the destructive effect of the central galaxy.
Differences between the curves become apparent outside the halo due to local
halo clustering. This effect originates predominantly from “backsplash" halos
that have pericenters with $r<R_{\rm vir}$ but apocenters of $r>R_{\rm vir}$,
meaning they lie within the splashback radius of the lens-mass host but spend
most of their time at $r>R_{\rm vir}$. For comparison, we plot the analytical
contribution of the two-halo term given by G18 (Eq. 1; dashed orange), where
we set $\delta_{\rm los}=1$. While G18 accurately reproduces the contribution
at $\ell>5\times R_{\rm vir}$ for most of the curves, it significantly under-
predicts the differential halo counts within $\ell\in[2-5]R_{\rm vir}$. Our
proposed modification, Eq. (3) (dashed black), captures the excess for
$\ell>2\,R_{\rm vir}$ to within $10\%$.

###  2.3 Counting within lens-centered cylinders

We quantify perturber counts by enumerating small halos that sit within
randomly-oriented lens-centered cylinders of length $\ell$ and radius
$\mathcal{R}$, where $\ell$ can extend over the length of the simulation box
and $\mathcal{R}$ will be set to a value close to an expected Einstein radius
($\sim 10$ kpc). As illustrated in Fig. 2, each cylinder is centered on a
target halo in the simulation box. The volume of the cylinder is increased by
lengthening $\ell\rightarrow\ell+\Delta\ell$ with the radius $\mathcal{R}$
fixed. For a discrete increment variation of $\Delta{\ell}$, both ends of the
cylinder are increased by $(\Delta\ell)/2$, which captures both structures
whose positions are found in the foreground and background of the target halo.
Doing so allows us to qualitatively compare counts as function of radius $r$
from the halo center, i.e., $r\simeq\ell/2$. For example, a lens-centered
cylinder length of $\ell=2R_{\rm vir}$ spans the full diameter of the dark
matter halo.

In what follows we compare counts from lens-centered cylinders to those of
average cylinders, which are configured like a lens-centered cylinder, but now
with their centers randomly placed within the simulation box. Note for a large
sample of average cylinders, the mean count per unit volume at any $\ell$ will
be equal to the average halo density (e.g. Sheth & Tormen, 2002). For both
cylinder types, we take into account the periodicity of the cosmological box
but never allow the full cylinder length to exceed the box length.

Fig. 3 provides a simple illustration of our counting prescription applied to
TNG100 for a single lens-mass halo of $M_{\rm vir}\simeq 10^{13}M_{\odot}$
(labeled as “lens-centered cylinder”) and for a randomly placed cylinder of
the same size inside the box (labeled as “average cylinder”).For illustrative
purposes, we have set the radius of the cylinder to a very large value
$\mathcal{R}=R_{\rm vir}\simeq 500\ \rm kpc$, much larger than what we will
use in our main analysis. The cylinder length is set to $\mathcal{\ell}=10$
Mpc. The filled circles show all subhalos with $M_{\rm halo}>10^{9}\
M_{\odot}$. Points and circles are colored based on their relative distance
from the host: within $\ell=2R_{\rm vir}$ ($r=R_{\rm vir}$) as cyan,
$\ell=[2-4]R_{\rm vir}$ ($r=[1-2]R_{\rm vir}$) as magenta, and black for
everything else out to a length $\ell=10$ Mpc. The two top plots show the
edge-on projections while the bottom plots are the face-on projections. For
the face-on projections, the sizes of the circles are scaled proportionally to
the mass of the halos in the cylinder (with the host halo removed for
clarity). In the bottom-left panel, we show only halos that are outside of
$R_{\rm vir}$ of the lens-host. Even excluding substructure, the overall count
is much higher than the random cylinder.

For this particular randomly-chosen "average" cylinder, we find 14 halos with
$M_{\rm halo}>10^{9}\ M_{\odot}$ within the 10 Mpc projection visualized. Note
that the cosmological average expected for this volume is $\approx 14.3$ when
using the Sheth & Tormen (2002) mass function. Counts are significantly higher
for the lens-centered cylinder. For comparison, the lens-centered cylinder
contains 108 halos of the same mass. Interestingly, correlated structure
counts around the lens-host that exist outside of $r=R_{\rm vir}$ but within
$r=2R_{\rm vir}$ is 28, which exceeds all of the counts within the 10 Mpc long
average cylinder. This shows that clustering in the vicinity of the lens will
boost signals non-trivially compared to what we would have estimated by
ignoring local clustering outside of the halo virial radius.

While Fig. 3 is useful to elucidate the point of a projected cylinder, the
radius shown, $\mathcal{R}=R_{\rm vir}$, is not relevant for lensing studies.
The remainder of the analysis hereafter imposes a fixed projected cylinder
radius of $\mathcal{R}=10\ \rm kpc$, which is a value comparable to size of
the lens-mass’ Einstein radius (typically $\lesssim$ 10 kpc). Note that in the
SLACS sample used by Vegetti et al. (2014) at a redshift $\langle z_{\rm
lens}\rangle\sim 0.2$, the median Einstein radius is $\sim 4.2$ kpc. We adopt
a slightly larger 10 kpc radius in what follows in order to improve counting
statistics (see Appendix A.1 for a discussion of counting variance). Our
primary results below are framed as relative counts per unit volume, such that
the precise radius of the cylinder factors out. We show in Appendix A.2 that
our results are insensitive (to within counting noise) to choices of cylinder
radii of $5$ kpc and even $2$ kpc.

##  3 Results

###  3.1 Average halo counts in projection

Our main results are presented in Fig. 4, where we plot the mean differential
count of halos per cylinder length, $\langle dN/d\ell\rangle$ using 100
cylinders randomly oriented around each lens-mass halo (solid lines).
Integrating $dN/d\ell$ over $\ell$ gives the cumulative count $N(<\ell)$
within a cylinder of total length $\ell$. We plot the differential count as a
function of cylinder length, $\ell$. Shown is the mean differential count
rather than the median because this allows us to compare directly to analytic
expectations for the average halo abundance. Each solid curve shows the rate
of counts for halo masses greater than a given value (indicated in the
legend), normalized by the rate of counts expected from the average background
from (Sheth & Tormen, 2002). Counts equal to the rate of counts from the
average background are shown as the horizontal dotted line with an amplitude
of $1$.555We have also tested with a large number average cylinders in the
simulations (as the example shown in Fig. 3) and confirmed that the average
background counts are consistent with analytical expectations of Sheth &
Tormen (2002). The vertical grey-dashed line separates between the two regimes
of interest: the substructure contribution $(\ell<2R_{\rm vir}$ and the local
structure $(\ell>2R_{\rm vir})$.

We see that for all simulations and mass cuts, the differential counts are
above the average counts out to $r\sim 20$ Mpc ($\ell\sim 40\,R_{\rm vir}$).
This is attributed to excess clustering in the vicinity of the massive host,
an effect often ignored in lensing studies (see e.g. Despali et al., 2018, who
assume average counts outside the virial radius), though some groups (e.g.
Gilman et al., 2018) have attempted to account for the effect (see below). For
perturbers more massive than $10^{9}\ M_{\odot}$, the rate of counts do not
reach the average background until $\ell\sim 75R_{\rm vir}\approx 40\ \rm Mpc$
for both TNG100 (magenta), TNG100DM (cyan), and FIREboxDM (red, which is
limited by only having three hosts). It is interesting also to compare TNG100
(magenta) to TNG100DM (cyan). We see that at small $\ell$ (corresponding to
the center of the host halo), the overall count is higher in the DMO run. This
comes about because of enhanced subhalo destruction from the central galaxy
potential (e.g. Graus et al., 2018). Differential counts in FIREboxDM are
consistent with those in TNG100DM for a lower mass threshold of $M_{\rm
halo}=10^{9}\,M_{\odot}$, though the FIREboxDM result is noisier owing to the
fact that there are only three target lens-mass halos in the volume. In
Appendix A.1, we further discuss the effect of sample variance in our
analysis.

Comparing the $10^{9}\ M_{\odot}$ (magenta/cyan), $10^{8}\ M_{\odot}$ (gray)
and $10^{7}\ M_{\odot}$ (black) lines, there is an indication that lower-mass
halos contribute more near the center of the lens-host ($\ell\sim 0.1R_{\rm
vir}$). This would be expected if subhalo radial distributions within the host
halo are more centrally concentrated at lower subhalo masses. Beyond the
virial radius $(\ell>2\,R_{\rm vir})$ the lower-mass halos found in FIREboxDM,
$10^{8}\ M_{\odot}$ (gray) and $10^{7}\ M_{\odot}$ (black), track the enhanced
counts seen at $10^{9}\ M_{\odot}$. For $\ell>5R_{\rm vir}$, the counts for
lower mass halos in FIREboxDM fall slightly below those seen in TNG. This
difference could be physical. For example, lower mass halos may be less
clustered around the lens host. However, the offset we see from $10^{9}\
M_{\odot}$ to $10^{8}\ M_{\odot}$ is much larger than would be expected
naively from the clustering bias change over this mass range:
$b(10^{9}M_{\odot})\simeq 0.64$ vs. $b(10^{8}M_{\odot})\simeq 0.63$ at $z=0.2$
(Sheth & Tormen, 1999). The difference could also arise from the lack of
large-scale power in the small volume of FIREboxDM or from simple sample
variance from having only three host halos. In Appendix A.1, we demonstrate
that these differences are consistent with sample variance; in order to
definitively confirm this, we would need a larger cosmological box with a
particle mass resolution comparable to FIREboxDM.

###  3.2 Analytic Model Comparison

Correlated structure outside of the virial radius of a massive target halo is
related to the "two-halo term" of the halo-matter correlation function (e.g.
Ma & Fry, 2000; Seljak, 2000; Smith et al., 2003). Gilman et al. (2018) (G18
hereafter) estimated this effect as

$\displaystyle\frac{d^{2}N}{dm\ dV}$ $\displaystyle=\delta_{\rm
los}\Big{[}1+\xi_{\rm 2halo}(r,M_{\rm vir},z)\Big{]}\frac{d^{2}N_{\rm ST}}{dm\
dV},\,,$ (1)

where $M_{\rm vir}$ is the mass of the host halo and $\delta_{\rm los}$ is an
overall scaling term that accounts for a systematic shift of the mean number
of halos predicted by the Sheth & Tormen (2002) mass function ($N_{\rm ST}$),
and

$\displaystyle\xi_{\rm 2halo}(r,M_{\rm vir},z)$ $\displaystyle=b_{\rm
ST}(M_{\rm vir},z)\xi_{\rm lin}(r,z)$ (2)

is the two-halo term that depends on the bias, $b_{\rm ST}$, around the lens
halo computed as in Sheth & Tormen (1999) and $\xi_{\rm lin}$ is the linear
matter-matter correlation function at a three-dimensional distance, $r$,
computed from the linear power spectrum at redshift $z$.

Figure 5:  Solid lines show the integrated mean count of halos more massive
than $10^{9}\ M_{\odot}$ within cylinders of radius $\mathcal{R}=10$ kpc and
of length $\ell$ centered on lens hosts from TNG100 (magenta) and TNG100DM
(cyan). The gray vertical line marks the cylinder length at the halo boundary
($\ell=2R_{\rm vir}$). The cumulative count of halos that would result without
any adjustment for correlated structure beyond the virial radius is shown by
the dotted lines. The predicted contribution of correlated structure from G18
(Eq. 3) is shown by the dashed lines. In implementing G18, we set $\delta_{\rm
los}=1$ for TNG100DM, as the convergence to average is about one-to-one.
However, for TNG100, we set $\delta_{\rm los}=0.879$ to remain consistent with
the average-count differences in TNG100DM. The G18 model captures much, but
not all, of the correlated structure; the extra component beyond the G18
prediction comes from virialized halos beyond $R_{\rm vir}$ but within the
splashback radius.

Eq. (1) (labeled “G18”) is plotted as the orange dashed curve in Fig. 4.
There, we designate $\delta_{\rm los}=1$ since the structure found in the
volumes of all three simulations tend to be well represented by the halo mass
function (refer to Fig. 1). Notice that in the region $\ell=[2-5]R_{\rm vir}$,
counts in the simulation are in excess of the G18 estimate. This excess likely
originates from subhalos with orbits that have apocenters beyond $R_{\rm
vir}$. This “backsplash" population can be substantial just outside of $R_{\rm
vir}$, as 50-80% of halos at $r\in[1-1.5]\,R_{\rm vir}$ were once subhalos
(e.g. Gill et al., 2005; Garrison-Kimmel et al., 2014) and therefore represent
a natural continuation of the subhalo population within $R_{\rm vir}$666A
backsplash excess is also hinted at pictorially in Fig. 3 (top panel) from the
magenta dots in the (edge-on) lens-centered projection. A physical boundary
for this virialized population of halos is the so-called “splashback" radius
of the host (More et al., 2015), where recently-accreted material reaches its
second apocenter (or the first apocenter after turn-around, where the turn-
around – or infall – radius is $R_{\rm infall}\approx 1.4\,R_{\rm sp}$). Our
sample of lens-mass halos at $z=0.2$ should have a median splashback radius of
$R_{\rm sp}\approx 1.5R_{\rm vir}$ (More et al., 2015). Subhalos outside of
$R_{\rm vir}$ but within $R_{\rm sp}$ ($\ell\approx 3R_{\rm vir}$), and
accompanying halos on first infall ($r<R_{\rm infall}\approx 2.1\,R_{\rm
vir}\leftrightarrow\ell\la 4.2\,R_{\rm vir}$) therefore provide a natural
explanation of the excess relative to G18 at $\ell=[2-5]R_{\rm vir}$. Notably,
the G18 model matches our simulation results for $\ell>5R_{\rm vir}$, similar
to the radius at which Garrison-Kimmel et al. (2014) find the backsplash
fraction is essentially zero.

To accommodate the excess clustering seen within $5R_{\rm vir}$, we modify Eq.
(1) as

$\displaystyle\frac{d^{2}N_{\rm CDM}}{dm\ dV}$ $\displaystyle=\delta_{\rm
los}\left[1+\xi_{\rm 2halo}(r,M_{\rm vir},z)+b_{\rm
sp}(r)\right]\,\frac{d^{2}N_{\rm ST}}{dm\ dV}\,,$ (3)

where

$\displaystyle b_{\rm sp}(r)$ $\displaystyle:=b_{e}\left(\frac{r}{5R_{\rm
vir}}\right)^{-s_{e}}\,.$ (4)

Here $b_{e}$ and $s_{e}$ are free parameters, and the term in parentheses
accounts for the excess clustering (compared to Eq. 1). Eq. (3) explicitly
separates the contribution of the cosmologically average LOS halos, an
enhancement from halo-halo clustering, and a further enhancement from
backsplash halos. The implied LOS count with this parametrization is shown as
the black-dashed curve in Fig. 4 with $b_{e}=0.1$ and $s_{e}=4$. This choice
of parameters captures our results to an accuracy of $10\%$ all the way out to
the edge of the simulation box. In particular, our parametrization is what we
would expect for a population of virialized (sub)halos that populate between
distances of $R_{\rm vir}$ and $5R_{\rm vir}$

###  3.3 Cumulative Counts

In order to estimate sample size of suitable lenses needed for testing
predictions, it is useful to know how many low-mass halos, on average, we
expect to see along the LOS to a host. The lower the expected count per lens,
the more lenses we will need to place meaningful constraints on the halo mass
function at low masses. The average count will depend, of course, on the
redshift of the target galaxy and lens relative to the observer (e.g. Despali
& Vegetti, 2017) but our results allow us to determine the average count
within the $\sim 100$ Mpc vicinity of the lens. Broadly speaking, the closer
the lens to the observer, the more important correlated structure will be. The
results that follow will be important for any lens where the substructure
contribution is significant compared to the total expected count of
perturbers. This is the case for roughly half of the lenses in the Despali &
Vegetti (2017) sample, for example.

Fig. 5 shows the mean count of halos along the projected lens-centered
cylinder as a function of cylinder length $\ell$, where the mean cumulative
counts, $\langle N(<\ell)\rangle$, is related to the average rate of counts by

$\displaystyle\left\langle
N(<\ell)\right\rangle=\int_{0}^{\ell}d\ell^{\prime}\
\left\langle\frac{dN}{d\ell^{\prime}}\right\rangle\,.$ (5)

Following the presentation of Fig. 4, for increasing $\ell$, the simulation
results fully realize the expected clustering contribution to the LOS halo
count once integrating the solid lines out to the cosmological boxes.

Solid lines in Fig. 5 show the integrated mean count of halos more massive
than $10^{9}\ M_{\odot}$ within cylinders of length $\ell$ centered on lens
hosts from TNG100 (magenta) and TNG100DM (cyan). The gray vertical line marks
the cylinder length at the halo boundary ($\ell=2R_{\rm vir}$). As we can see
by comparing TNG100 and TNG100DM at $\ell<2R_{\rm vir}$, subhalo counts are
reduced by $\approx 40$ in the full-physics run compared to the DMO run. This
result is consistent with the findings of Despali & Vegetti (2017) and Graus
et al. (2018). A second takeaway from Fig. 5 is that the average count of
halos is much less than unity out to the edge of the box. This means that most
of the lens-centered LOS cylinder do not contain halos larger than $10^{9}\
M_{\odot}$within within 100 Mpc in projections of radius $\mathcal{R}=10$kpc.
We find that for TNG100DM, $87.5\%$ of projections contain no halos at this
mass limit. For TNG100, the fraction of empty projections rises to $\approx
92.5\%$.

Figure 6:  Similar to Fig. 5, but now for FIREboxDM halos with $M_{\rm
halo}>10^{8}\ M_{\odot}$ (top) and $M_{\rm halo}>10^{7}\ M_{\odot}$ (bottom).
Both curves (top and bottom panels) incorporating the G18 model use
$\delta_{\rm los}=1$. The typical LOS passing through a lensing host
encounters significantly more perturbers above a given mass as the perturber
mass threshold is lowered.

The dashed and dotted lines in Fig. 5 compares our simulation results to
alternative ways of estimating cumulative counts beyond the host halo virial
radius. The first assumes that halo counts are equal to the universal average
(estimated via Sheth & Tormen, 2002) for all radii beyond the virial radius of
the lens (dotted lines) and the second uses the estimate from G18, which
models local clustering as in Eq. 1. We see that if we assume the average
background is achieved for radii beyond the virial radius, the cumulative
count at 100 Mpc is underpredicted by $\sim 20\%$ in TNG100 and $\sim 15\%$ in
TNG100DM. Differences between the simulations and the G18 estimate are not as
large, with $\approx 5\%$ offsets in TNG100DM and TNG100 at 100 Mpc. Note that
when using the G18 formula, we set $\delta_{\rm los}=1$ for TNG100DM and
$\delta_{\rm los}=0.879$ for TNG100. The latter value is below unity because
the average differential count is slightly below the Sheth & Tormen (2002)
estimate at large $\ell$ (see the magenta vs. dotted lines at > 40 Mpc in Fig.
4). This factor is also introduced for the average background interpolation.
We provide a more thorough discussion on how the contribution of the
clustering component to quantifying the LOS structure of TNG100 in later in
this section, while TNG100DM results are discussed in Appendix B.

Fig. 6 displays the mean cumulative counts for FIREboxDM subhalos of $M_{\rm
halo}>10^{8}\ M_{\odot}$ and $M_{\rm halo}>10^{7}\ M_{\odot}$ in the top and
bottom panel, respectively. The likelihood of finding a halo in a single
projected cylinder increases substantially at these lower masses compared to
the $10^{9}\ M_{\odot}$. Specifically, we expect to see, on average, more than
one small halo per LOS for $10^{8}\ M_{\odot}$ and more than 15 for $10^{7}\
M_{\odot}$.One caveat here is that these simulations do not include the
destructive effects of a central galaxy. We would expect substructure to be
depleted within $\ell\approx 1$ Mpc by approximately $40\%$ if the destruction
mirrors that seen in Fig. 5 for $M_{\rm halo}>10^{9}\ M_{\odot}$
halos.777Recently, Kelley et al. (2019) presented high-resolution zoom
simulations for MW-mass DMO halos while accounting for the central galaxy and
found that the depletion of substructure is about roughly the same factor for
halos down to $10^{7}\ M_{\odot}$ at $z=0$. While we are drawing possible
conclusions from MW-mass halos, this scaling could translate to our lens-mass
halos.

The line styles in Fig. 6 mimic those in Fig. 5, with solid lines representing
the full simulation results. The line is extrapolated beyond the edge of
FIREboxDM (gray solid line) by assuming it follows the G18 estimate starting
at $L_{\rm box}^{\rm FIRE}$ with $\delta_{\rm los}=1$ for both mass cuts. We
have also modeled the entire solid line using Eq. (3) instead of the
simulation results to test whether the relatively small box of FIREboxDM
suppresses large-scale modes, thereby affecting the number of halos found from
$\ell=5R_{\rm vir}$ out to $L_{\rm box}^{\rm FIRE}$. Doing so, we find very
little difference in the final results. The assumption of average background
to the simulation results out to $L_{\rm box}^{\rm TNG}$ results differ to
about $15\%$ for $10^{8}\ M_{\odot}$ while $10\%$ for $10^{7}\ M_{\odot}$.

It is important to quantify the contribution of correlated clustering to
counts for $\ell>2\,R_{\rm vir}$, i.e., when subtracting off the contribution
of subhalos to $N(<\ell)$, what is the fractional contribution of the
population of clustered halos compared to the average halo population at
$\ell>2\,R_{\rm vir}$? This is explicitly address in Fig. 7, which shows the
total average cumulative count of $10^{9}\ M_{\odot}$ halos in TNG100 for
cylinder lengths ranging from the edge of the virial radius to the edge of the
box. We plot the average count $\langle N(<\ell)\rangle$ in units of the
average cumulative count from subhalos $\langle N_{\rm
sub}\rangle\equiv\langle N(\ell<2R_{\rm vir})\rangle$, where $\langle N_{\rm
sub}\rangle=0.072$. The full count from the simulation (labeled as “total”) is
shown by the magenta curve while cylinders that are assumed to have the
average background outside of the virial radius (labeled “average”) are
plotted as the light gray curve. To quantify the excess clustering associated
beyond substructure, we take the difference between the magenta and gray
curves. This results in the black curve (“clustering”). The excess clustered
contribution asymptotes to $\sim 1.35\langle N_{\rm sub}\rangle$ at
$\ell\approx 70\ \mathrm{Mpc}\ ({\sim}75\,R_{\rm vir})$. This means that local
clustering boosts the expected signal by $\sim 35\%$ compared to what we would
expect from subhalos alone. We show similar results for lower mass halos in
other simulations in Appendix B. Broadly speaking, the boost from local
clustering is smaller relatively in DMO simulations that do not have enhanced
subhalo destruction from the central galaxy. Note that the clustering
contribution is larger than the average contribution out to $\ell\approx
65\,{\rm Mpc}$, or equivalently, ${\sim}130\,R_{\rm vir}$.

Figure 7: — Boost in counts from local clustering: The mean count of LOS halos
more massive than $10^{9}\ M_{\odot}$ within lens-centerd cylinders from
TNG100 relative to the mean count within the virial radius $\langle N_{\rm
sub}\rangle\equiv\langle N(\ell<2\,R_{\rm vir})\rangle$ is presented here. The
magenta line shows the mean total count from the simulation (labeled “total”),
while the gray line shows counts in cylinders that assume the density of halos
matches the mean background (dubbed “average”) for $\ell>2\,R_{\rm vir}$. The
black line shows the difference between the “total" and “average"
contributions; this is the component that is attributable to local halo
clustering (dubbed “clustering”). We see that the local clustering effect
provides a boost of $\sim 35\%$ compared to the subhalo count alone and that
this contribution dominates the “average" contribution out to ${\sim}65\,{\rm
Mpc}\approx 130\,R_{\rm vir}$.

###  3.4 Structure along principal axes

CDM halos have significant triaxiality (e.g. Frenk et al., 1988; Dubinski &
Carlberg, 1991; Warren et al., 1992; Cole & Lacey, 1996; Jing & Suto, 2002;
Bailin & Steinmetz, 2005; Kasun & Evrard, 2005; Paz et al., 2006; Allgood et
al., 2006; Bett et al., 2007; Muñoz-Cuartas et al., 2011; Despali et al.,
2014; Vega-Ferrero et al., 2017; Lau et al., 2020), which could also mean that
subhalos are found preferentially along the host’s major (densest) axis (e.g.
Zentner et al., 2005). This effect can qualitatively impact our analysis of
lens-centered projections, especially for substructure lensing. Furthermore,
it is likely that galaxy lenses are biased to be oriented with the LOS
coinciding with the host halo’s major axis (Mandelbaum et al., 2009; Osato et
al., 2018), as this configuration produces a larger surface mass density for a
fixed overall mass distribution. In order to explore the potential magnitude
of this effect, we calculate dark matter halo shapes from the lens targets
using the shape inertia tensor as discussed in Allgood et al. (2006). We
include all dark matter particles within a shell between $10-20\%$ of $R_{\rm
vir}$ as a conservative approach for our sample of halos. Using a shell rather
than the enclosed mass minimizes the influence of particles with radii smaller
than the numerical convergence scale. The resulting eigenvalues of the shape
tensor are proportional to the square root of the principal axes of the dark
matter distribution. We then re-do the analysis of Section 3, now aligning the
lens-centered projections along each principal axis of the lens-mass halo.

Figure 8: — Projections along the principal axes: Similar presentation as Fig.
4, but now filled solid lines show the counts for TNG100 halos $M_{\rm
halo}>10^{9}\ M_{\odot}$ along the major, intermediate, and minor axis in the
top, center, and bottom panels, respectively. The transparent lines give the
mean counts of the substructure with local clustering counts of TNG100
presented previously in Fig. 5. Sight-lines oriented along the major (minor)
axis of the halo result in a non-trivial increase (decrease) in the number of
perturbers encountered along the sight line.

Fig. 8 plots the results for the average (mean) counts along each principal
axis for halos more massive than $10^{9}\ M_{\odot}$ in TNG100 (for comparison
with TNG100DM, see Fig. 10). The top, middle, and bottom panels depict the
average counts along the major axis, intermediate axis, and minor axis,
respectively, with cylinders of radius $\mathcal{R}=10\ \rm kpc$. For
comparison, the faded solid line in each panel shows the mean counts presented
previously in Fig. 5. We see that the major axis sight-line results in
measurably more halos, on average, than do other orientations. The boost along
the major axis is $\sim 30\%$ compared to the random LOS, with essentially the
entire contribution coming from subhalos (and backsplash halos). Along the
minor axis, on the other hand, average counts are significantly reduced.
Projections along the intermediate-axis are comparable to the average counts
shown in Fig. 4. It is clear that lens-centered projections along the major or
minor axis can non-trivially boost or decrease the lensing signals by a factor
of $\sim 2$. This factor is also acquired for analog simulation neglecting the
presence of the central galaxy and baryons (see Fig. 10).

##  4 Summary

Using the set of TNG100 and DMO FIREboxDM cosmological simulations, we
quantify the effect of local clustering on gravitational lensing searches for
low-mass dark matter halos. We specifically focus on lens-mass halos of mass
$M_{\rm vir}\simeq 10^{13}\ M_{\odot}$ at $z=0.2$ as prime targets for future
lensing surveys and explore counts of halos down to $M_{\rm halo}=10^{7-9}\
M_{\odot}$.

Our primary result is that local clustering can boost the expected LOS
perturber halo counts significantly compared subhalos alone. The signal
exceeds that expected for an average background projection to distances in
excess of $\pm 10$ Mpc from the lens host (Fig. 4), with a significant excess
within $2-5\,R_{\rm vir}$. We provide an analytic expression for this
contribution (Eq. 3), which we hope will be useful in full lensing
interpretation studies.

Using full-physics TNG100 simulations (which resolve halos down to $10^{9}\
M_{\odot}$), we find that the central galaxy in lens-mass hosts depletes
subhalos by $\sim 70\%$ compared to dark-matter-only simulations (see Fig. 5).
This result agrees with previous studies (Despali & Vegetti, 2017; Graus et
al., 2018). From TNG100, the excess local clustering outside of the virial
radius gives an expected count that is $\sim 35\%$ higher than would be
expected from subhalos alone (Fig. 7).

Local contributions to perturber counts are also affected by halo orientation.
The above results assume a random lens-orientation with respect to the
observer, but if there is a bias for lenses to be oriented along the principal
axis (e.g. Dietrich et al., 2014; Groener & Goldberg, 2014), then the expected
local count may be enhanced. Our initial exploration of this issue indicates
that local projected counts are $\sim 50\%$ higher when the target halo is
oriented along the major axis compared to a random orientation (Fig. 8). This
result, in turn, has implications for derived constraints on the mass spectrum
of perturbers and accompanying constraints on dark matter particle properties.

The above results will be particularly important for low-redshift lenses
($z_{l}\sim 0.2$), such as those in the SLACS sample (Vegetti et al., 2014).
For such lenses, Despali et al. (2018) find that subhalos should contribute
$\sim 30-50\%$ of the total perturber signal relative to LOS halos that
neglect local clustering. With clustering included, the local (subhalo +
clustering) contribution may well be comparable to $\sim 67.5\%$ of the total
LOS contribution that neglect clustering (or $\sim 40\%$ of the total
contribution) for some lenses, especially those with lower-redshift sources
($z_{s}\sim 0.6$). Taking into account any biases in lens-halo orientation is
also important for these lenses, as an additional $\sim 50\%$ boost in counts
from local clustering could significantly affect interpretations.

## Acknowledgements

This article was worked to completion during the COVID-19 lock-down and would
not have been possible without the labors of our essential workers. We thank
Ran Li for helpful comments in improving the early version of this article. We
are thankful to Quinn Minor and Manoj Kaplinghat for helpful discussions. The
authors thank the Illustris collaboration for facilitating the IllustrisTNG
simulations for public access and the FIRE collaboration for the use of the
DMO FIREbox simulation for our analysis. AL and JSB was supported by the
National Science Foundation (NSF) grant AST-1910965. AL was partially
supported by NASA grant 80NSSSC20K1469. MBK acknowledges support from NSF
CAREER award AST-1752913, NSF grant AST-1910346, NASA grant NNX17AG29G, and
HST-AR-15006, HST-GO-14191, HST-GO-15658, HST-GO-15901, and HST-GO-15902 from
the Space Telescope Science Institute, which is operated by AURA, Inc., under
NASA contract NAS5-26555. RF acknowledges financial support from the Swiss
National Science Foundation (grant no 157591 and 194814). The analysis in this
manuscript made extensive use of the python packages COLOSSUS (Diemer, 2018),
NumPy (van der Walt et al., 2011), SciPy (Oliphant, 2007), and Matplotlib
(Hunter, 2007); We are thankful to the developers of these tools. This
research has made all intensive use of NASA’s Astrophysics Data System
(http://ui.adsabs.harvard.edu/) and the arXiv eprint service
(http://arxiv.org).

## Data Availability

The data supporting the plots within this article are available on reasonable
request to the corresponding author. A public version of the Arepo code is
available at https://arepo-code.org/about-arepo. The IllustrisTNG data,
including simulation snapshots and derived data products are publicly
available at https://www.tng-project.org/. A public version of the GIZMO code
is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.
Additional data from the FIRE project, including simulation snapshots, initial
conditions, and derived data products, are available at
https://fire.northwestern.edu/data/.

## References

  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
  * Bailin & Steinmetz (2005) Bailin J., Steinmetz M., 2005, ApJ, 627, 647
  * Banik et al. (2019) Banik N., Bovy J., Bertone G., Erkal D., de Boer T. J. L., 2019, arXiv e-prints, p. arXiv:1911.02663
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
  * Bertone & Tait (2018) Bertone G., Tait T. M. P., 2018, Nature, 562, 51
  * Bett et al. (2007) Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
  * Bode et al. (2001) Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93
  * Bolton et al. (2008) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 2008, ApJ, 682, 964
  * Bond et al. (1996) Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603
  * Bozek et al. (2016) Bozek B., Boylan-Kolchin M., Horiuchi S., Garrison-Kimmel S., Abazajian K., Bullock J. S., 2016, MNRAS, 459, 1489
  * Bradač et al. (2002) Bradač M., Schneider P., Steinmetz M., Lombardi M., King L. J., Porcas R., 2002, A&A, 388, 373
  * Brooks & Zolotov (2014) Brooks A. M., Zolotov A., 2014, ApJ, 786, 87
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Bullock et al. (2000) Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517
  * Carlberg (2009) Carlberg R. G., 2009, ApJ, 705, L223
  * Cole & Lacey (1996) Cole S., Lacey C., 1996, MNRAS, 281, 716
  * Collett (2015) Collett T. E., 2015, ApJ, 811, 20
  * D’Aloisio et al. (2014) D’Aloisio A., Natarajan P., Shapiro P. R., 2014, MNRAS, 445, 3581
  * Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Despali & Vegetti (2017) Despali G., Vegetti S., 2017, MNRAS, 469, 1997
  * Despali et al. (2014) Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208
  * Despali et al. (2018) Despali G., Vegetti S., White S. D. M., Giocoli C., van den Bosch F. C., 2018, MNRAS, 475, 5424
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Dietrich et al. (2014) Dietrich J. P., et al., 2014, MNRAS, 443, 1713
  * Dolag et al. (2016) Dolag K., Komatsu E., Sunyaev R., 2016, MNRAS, 463, 1797
  * Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
  * Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948
  * Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
  * Feldmann & Spolyar (2015) Feldmann R., Spolyar D., 2015, MNRAS, 446, 1000
  * Frenk et al. (1988) Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
  * Garrison-Kimmel et al. (2014) Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K., 2014, MNRAS, 438, 2578
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Geller & Huchra (1989) Geller M. J., Huchra J. P., 1989, Science, 246, 897
  * Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  * Gill et al. (2005) Gill S. P. D., Knebe A., Gibson B. K., 2005, MNRAS, 356, 1327
  * Gilman et al. (2017) Gilman D., Agnello A., Treu T., Keeton C. R., Nierenberg A. M., 2017, MNRAS, 467, 3970
  * Gilman et al. (2018) Gilman D., Birrer S., Treu T., Keeton C. R., Nierenberg A., 2018, MNRAS, 481, 819
  * Gilman et al. (2019) Gilman D., Birrer S., Treu T., Nierenberg A., Benson A., 2019, MNRAS, 487, 5721
  * Graus et al. (2018) Graus A. S., Bullock J. S., Boylan-Kolchin M., Nierenberg A. M., 2018, Monthly Notices of the Royal Astronomical Society, 480, 1322
  * Groener & Goldberg (2014) Groener A. M., Goldberg D. M., 2014, ApJ, 795, 153
  * He et al. (2018) He Q., Li R., Lim S., Frenk C. S., Cole S., Peng E. W., Wang Q., 2018, MNRAS, 480, 5084
  * He et al. (2020) He Q., et al., 2020, arXiv e-prints, p. arXiv:2010.13221
  * Hezaveh et al. (2016) Hezaveh Y. D., et al., 2016, ApJ, 823, 37
  * Horiuchi et al. (2016) Horiuchi S., Bozek B., Abazajian K. N., Boylan-Kolchin M., Bullock J. S., Garrison-Kimmel S., Onorbe J., 2016, MNRAS, 456, 4346
  * Hsueh et al. (2017) Hsueh J. W., et al., 2017, MNRAS, 469, 3713
  * Hsueh et al. (2018) Hsueh J.-W., Despali G., Vegetti S., Xu D. a., Fassnacht C. D., Metcalf R. B., 2018, MNRAS, 475, 2438
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Ibata et al. (2002) Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915
  * Inoue (2016) Inoue K. T., 2016, MNRAS, 461, 164
  * Inoue & Takahashi (2012) Inoue K. T., Takahashi R., 2012, MNRAS, 426, 2978
  * Inoue et al. (2015) Inoue K. T., Takahashi R., Takahashi T., Ishiyama T., 2015, MNRAS, 448, 2704
  * Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
  * Kasun & Evrard (2005) Kasun S. F., Evrard A. E., 2005, ApJ, 629, 781
  * Kelley et al. (2019) Kelley T., Bullock J. S., Garrison-Kimmel S., Boylan-Kolchin M., Pawlowski M. S., Graus A. S., 2019, MNRAS, 487, 4409
  * Khandai et al. (2015) Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS, 450, 1349
  * Koopmans (2005) Koopmans L. V. E., 2005, MNRAS, 363, 1136
  * Lau et al. (2020) Lau E. T., Hearin A. P., Nagai D., Cappelluti N., 2020, MNRAS,
  * Li et al. (2016) Li R., Frenk C. S., Cole S., Gao L., Bose S., Hellwing W. A., 2016, MNRAS, 460, 363
  * Li et al. (2017) Li R., Frenk C. S., Cole S., Wang Q., Gao L., 2017, MNRAS, 468, 1426
  * Ma & Fry (2000) Ma C.-P., Fry J. N., 2000, ApJ, 538, L107
  * MacLeod et al. (2013) MacLeod C. L., Jones R., Agol E., Kochanek C. S., 2013, ApJ, 773, 35
  * Mandelbaum et al. (2009) Mandelbaum R., van de Ven G., Keeton C. R., 2009, MNRAS, 398, 635
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Meneghetti et al. (2020) Meneghetti M., et al., 2020, Science, 369, 1347
  * Metcalf (2005) Metcalf R. B., 2005, ApJ, 629, 673
  * Metcalf & Amara (2012) Metcalf R. B., Amara A., 2012, MNRAS, 419, 3414
  * More et al. (2015) More S., Diemer B., Kravtsov A. V., 2015, ApJ, 810, 36
  * Muñoz-Cuartas et al. (2011) Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
  * Nadler et al. (2021) Nadler E. O., Birrer S., Gilman D., Wechsler R. H., Du X., Benson A., Nierenberg A. M., Treu T., 2021, arXiv e-prints, p. arXiv:2101.07810
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nierenberg et al. (2017) Nierenberg A. M., et al., 2017, MNRAS, 471, 2224
  * Oliphant (2007) Oliphant T. E., 2007, Computing in Science and Engineering, 9, 10
  * Osato et al. (2018) Osato K., Nishimichi T., Oguri M., Takada M., Okumura T., 2018, MNRAS, 477, 2141
  * Paz et al. (2006) Paz D. J., Lambas D. G., Padilla N., Merchán M., 2006, MNRAS, 366, 1503
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 475, 648
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Richings et al. (2021) Richings J., Frenk C., Jenkins A., Robertson A., Schaller M., 2021, MNRAS, 501, 4657
  * Sánchez et al. (2006) Sánchez A. G., Baugh C. M., Percival W. J., Peacock J. A., Padilla N. D., Cole S., Frenk C. S., Norberg P., 2006, MNRAS, 366, 189
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schneider et al. (2013) Schneider A., Smith R. E., Reed D., 2013, MNRAS, 433, 1573
  * Seljak (2000) Seljak U., 2000, MNRAS, 318, 203
  * Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
  * Sheth & Tormen (2002) Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
  * Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Tegmark et al. (2004) Tegmark M., et al., 2004, ApJ, 606, 702
  * Vega-Ferrero et al. (2017) Vega-Ferrero J., Yepes G., Gottlöber S., 2017, MNRAS, 467, 3226
  * Vegetti et al. (2010) Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969
  * Vegetti et al. (2014) Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014, MNRAS, 442, 2017
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Warren et al. (1992) Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weiner et al. (2020) Weiner C., Serjeant S., Sedgwick C., 2020, Research Notes of the American Astronomical Society, 4, 190
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23
  * Xu et al. (2009) Xu D. D., et al., 2009, MNRAS, 398, 1235
  * Xu et al. (2012) Xu D. D., Mao S., Cooper A. P., Gao L., Frenk C. S., Angulo R. E., Helly J., 2012, MNRAS, 421, 2553
  * Xu et al. (2015) Xu D., Sluse D., Gao L., Wang J., Frenk C., Mao S., Schneider P., Springel V., 2015, MNRAS, 447, 3189
  * Yoon et al. (2011) Yoon J. H., Johnston K. V., Hogg D. W., 2011, ApJ, 731, 58
  * Zentner et al. (2005) Zentner A. R., Kravtsov A. V., Gnedin O. Y., Klypin A. A., 2005, ApJ, 629, 219
  * Zhu et al. (2016) Zhu Q., Marinacci F., Maji M., Li Y., Springel V., Hernquist L., 2016, MNRAS, 458, 1559
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22

##  Appendix A Counting statistics

###  A.1 Sampling variance

In Section 3.1, we discussed the possibility that the differences in local
clustering signal between TNG100DM and FIREboxDM could arise from sampling
variance – there are only three lens-mass hosts in FIREboxDM – rather than
differences in large-scale structure (compare the red and cyan curve in Fig.
4). In order to explore this more fully, Fig. 9 shows $\langle
dN/d\ell\rangle$ for TNG100DM (cyan) and FIREboxDM (red). This figure mirrors
Fig. 4 except now we are using a cylinder radius with the size of the lens-
mass virial radius ($\mathcal{R}=R_{\rm vir}\simeq 500$ kpc) as opposed to the
typical Einstein radius ($\mathcal{R}=10$ kpc) in order to improve counting
statistics. The solid lines depict the mean of all projections while the
dashed (shown only for TNG100DM) is the median of all projections. The cyan
band encloses the $\pm 1\sigma$ of all of the individual curves around median.

Figure 9:  Analogous to Fig. 4, now with results of TNG100DM (cyan) and
FIREboxDM (red) for $M_{\rm halo}>$ $10^{9}\ M_{\odot}$ for cylinders with a
radius size comparable to the lens mass halo ($\mathcal{R}=R_{\rm vir}$) as
opposed to the lens halo’s Einstein radius ($\mathcal{R}=10\,{\rm kpc}$). The
solid curves depict the mean differential count while the dashed cyan is the
_m_ edian differential count for TNG100DM. The cyan band encloses the $\pm
1\sigma$ region of the curves. The median counts fall below the mean counts,
indicating that the mean counts are strongly affected by rare, high-count
orientations.

For $\ell<1$ Mpc, the median and mean curve for TNG100DM are consistent with
one another, but for $\ell>1$ Mpc we see that the average count (which is what
we plot in Fig. 4) sits well above the median. This is indicative of a highly
non-Gaussian distribution, with a tail skewed towards large fluctuations, as
expected for non-linear dark matter structure. Specifically, rare (high-count)
events drive the average higher than the median. If this is the case, we
require a large number of realizations in order to sample enough of the
distribution to capture the true average. As further evidence that sample
variance is the cause of the differences in local clustering signal between
TNG100DM and FIREboxDM, the median line in Fig. 9 falls below the average
background (black dashed) at $\ell\approx 10$ Mpc. This is due to a large
number of null counts.

Though the average count from FIREboxDM falls below the average line from
TNG100DM, we see that it falls within the $1\sigma$ band about the median.
Given that we only have three host halos in our sample, this would not be
unexpected even if the halos were sampling the same large-scale structure
field. Naively, there is a $\sim 30\%$ chance that three randomly drawn
distributions will lie within $1\sigma$ of the median. Given that we only have
three host halos, we conclude that the observed result is consistent with
small-number statistics.

Figure 10:  Analogous to Fig. 4, now for cylinders with a radius size
$\mathcal{R}=5$ kpc and a much smaller annulus $\mathcal{R}=2$ kpc. While the
curves are much noisier as the radius size becomes smaller (owing to the
increased fraction of null sight lines), the we are able to recover the same
trend for the differential counts seen for a projected radius of
$\mathcal{R}=10$ kpc.

###  A.2 Choice of cylinder radius

The conclusions made in the main text are based on the average number density
of subhalos within a project cylinder radius of 10 kpc, which is comparable to
the typical Einstein ring of our $10^{13}\ M_{\odot}$ lens-mass halos. The
choice of $\mathcal{R}=10\rm\,kpc$ could bias our results towards higher
substructure, and potentially, LOS halo counts, as massive galaxy lenses
typically have Einstein ring radii of $\sim 5$ kpc (Bolton et al., 2008).

In Fig. 10, we demonstrate how the selection of smaller lens radii does not
impact the trends discussed in our main analysis. Both the top and bottom plot
mirrors exactly Fig. 4, now with the top and bottom plot showing results for a
projected radius of 5 and 2 kpc, respectively. We find that we are able to
mostly recover the same trend seen for the 10 kpc projected radius from Fig. 4
for both $\mathcal{R}=$ 5 and 2 kpc. In the 5 kpc radius projection, the
curves are more noisier than the 10 kpc owing to the increasing fraction of
projections with zero halos. We would argue that the $\mathcal{R}=5$ kpc size
projections have comparable trends to the $\mathcal{R}=10$ kpc case. Out to
$\ell=5R_{\rm vir}$, most of the mean differentials drives faster down to the
average compared to the 10 kpc case, but we suspect this is owed to low number
statistics along with the sampling variance for TNG100DM (see discussion in
Appendix A.1). This becomes more apparent for the smaller projected radius of
2 kpc, as the fraction of projections of zero halos increases appreciably.
Note that the $M_{\rm halo}>$$10^{9}\ M_{\odot}$ result is impacted greatly by
sampling variance the low number of halos along the LOS.

##  Appendix B Supplementary Discussion with Dark Matter Only Physics

This appendix presents the DMO results from TNG100DM and for FIREboxDMand
provides additional elaboration of the results in the main text. Section B.1
picks up from later discussion of Section 3.3 and Section B.2 picks up from
Section 3.4.

Figure 11:  Like Fig. 7, now showing only the results for TNG100DM and
FIREboxDM. In the top panel, the black and gray curve are computed based off
of the total counts from TNG100DM (cyan curve). As a comparison, the blue
curve shows the total counts from FIREboxDM out to $L^{\rm FIRE}_{\rm box}$.
The middle and bottom panel are presented similarly with their curves
extrapolated out to the $L^{\rm TNG}_{\rm box}$ using Eq. (3).

###  B.1 Clustering contribution to the line-of-sight population

Fig. 11 presents the contribution of the simulated clustering to the LOS halo
population in TNG100DM and FIREboxDM. Namely, the top, middle, and bottom
panels plots the subhalo populations based on the lower-mass cuts of $10^{9}\
M_{\odot}$, $10^{8}\ M_{\odot}$, and $10^{7}\ M_{\odot}$, respectively. The
colored curves, labeled “total”, are the actual results of the simulation out
to the box using our method discussed in the main text. The average background
expectation for halos above the low-mass cuts in each panel is depicted by the
grey curves. The clustering contribution from the simulations, labeled
“clustering”, are plotted as the black curves, which is the difference between
the “total” curve and “average” curves. Like before, all results are
normalized by the $\langle N_{\rm sub}\rangle$ quantified from out lens-mass
sample. Since the actual FIREboxDM results extend out to $L_{\rm box}^{\rm
FIRE}$, we extrapolate the curves out to $L_{\rm box}^{\rm TNG}$using Eq. (1)
(using Eq. 3 is also equally adequate).

Figure 12:  Like Fig. 8, but now showing the TNG100DM counts for halos $M_{\rm
halo}>10^{9}\ M_{\odot}$ along the major, intermediate, and minor axis in the
top, center, and bottom panels, respectively. Plotted for comparison, the
transparent lines are the mean counts of the substructure with local
clustering counts for TNG100DM presented previously in Fig. 5.

Starting with the counts $10^{9}\ M_{\odot}$ (top panel), clustering acts as
the most contributing component to the LOS halo population until the average
background takes over at $\ell\approx 75R_{\rm vir}$ (or $r\approx 37.5\ \rm
Mpc$) for TNG100DM (cyan curve). As a comparison check, we plotted the
FIREboxDM results of $10^{9}\ M_{\odot}$ (thin blue curve) and find minimal
difference in results. Note the clustering component is determined from
TNG100DM and not FIREboxDM. The clustering component in TNG100DM boosts about
$1.5\langle N_{\rm sub}\rangle$ out until the average background takes over at
$\ell\approx 70\ \rm Mpc$ (or $r\approx 35$ Mpc). This is less significant
than the TNG100 in the main text, the clustering component strongly boosts the
number of LOS halos to about $1.7\langle N_{\rm sub}\rangle$ until the average
background takes over.

As we go further down to the low-mass cuts of the subhalo populations, the
contribution from correlated clustering becomes much weaker for decrease mass.
For $M_{\rm halo}>10^{8}\ \rm M_{\odot}$ (middle panel), the clustering
component boosts the number of LOS halos to only about $\sim 1.2\langle N_{\rm
sub}\rangle$ until the average background takes over at about $\ell\approx 60\
\rm Mpc$ (or $r\approx 30$ Mpc). The clustering component for $M_{\rm
halo}>10^{7}\ \rm M_{\odot}$ (bottom panel) becomes weaker by only boosting
the number of LOS halos to only about $\sim 1.1\langle N_{\rm sub}\rangle$
until the average background takes over at $r\approx 30$ Mpc). In order to
conclude with more robust predictions, a larger sample lens-mass halos in
comparable cosmological environments is needed to reduced the uncertainty
possibly found in $\langle N_{\rm sub}\rangle$ for these low-mass halos.

###  B.2 Structure along principal axes

Fig. 12 plots the mean subhalo counts as seen from lens-centered projections
along the principal axis of the TNG100DM lens-mass halos. The axes were
computed using the method detailed in Section 3.4. Shown are only the results
from TNG100DM since these cosmological volumes have enough lens-centered mass
halos to provide enough statistics since each halo has only three principal
axis to orientate on. Moreover, FIREboxDM, while useful for simulating
$10^{7-8}\ M_{\odot}$ halos in cosmological environment, has only three lens-
centered hosts, which will not provide adequate statistics to present. Though,
we find similar trends. The top, middle, and bottom panels depict the average
counts along the major-, intermediate-, and minor-axis, respectively. For
comparison, the faded solid line in each panel plots the mean counts presented
previously in Fig. 5. We again see that projections along the major-axis of
DMO lens-mass halos, subhalos tend to populate on the densest principal axis.
This effect is more dramatic for the substructure compared to the TNG100 shown
previously, owing to the lack of a central galaxy.

Projections the size of the Einstein radius along the major-axis find around
$75\%$ of the projections to contain no substructure out to the size of the
halo. Along the minor-axis, this is around $\sim 90\%$. Additionally, the
substructure counts appear to be almost comparable to the mean from Fig. 4,
though slightly less. Done on the intermediate-axis, this results in counts
somewhat comparable to the mean counts, though we see ${\sim}5\%$ boost in the
LOS component.

##  Appendix C Comparison between mass functions of lens-mass sample

As mentioned in Section 2, the SUBFIND halo finder was ran for the TNG suites
while ROCKSTAR was applied to FIREboxDM. Fig. 1 provides confidence in the
resulting halo catalogs used in the main text, as both catalogs show robust
agreement with the analytical prediction. However, halo finders vary in
routines for quantifying the masses for subhalos and can potentially produce
different mass functions when ran to the same simulations. It would be worth
comparing our resulting halo finders with one another based on the resolved
substructure population found for our lens-target halos at $z=0.2$.

Fig. 13 shows the subhalo mass functions for the $10^{13}\ M_{\odot}$ target-
lens systems in TNG100DM (purple) while the three from FIREboxDM (black
curves). The purple band for TNG100DM encloses the $90\%$ dispersion of all of
the target-lens subhalo mass functions. Significant disagreement could arise
based on the assignment of particles to subhalos from halo finders, (Graus et
al., 2018). To clarify this point, Fig. 14 provides a different view of
subhalo sample, where the subhalo $V_{\rm max}$ function for the same lens-
target systems.

Figure 13:  The subhalo mass functions for the lens-target halos in our main
sample compared between the TNG100DM (purple curve) and FIREboxDM (black
curves). FIREboxDM follows the median TNG100DM curves sufficiently for the
subhalo masses ranging from $\sim 10^{9-10}\ M_{\odot}$ while the three lens-
target in FIREboxDM hosts several more massive $10^{11}\ M_{\odot}$. Though,
the curves are still mostly found enclosed in the $90\%$ dispersion. The
agreement depicted from the subhalo mass functions provides further
justification for the choice of our ROCKSTAR parameters detailed in Section 2.
Figure 14:  The subhalo $V_{\rm max}$ functions for the lens-target halos in
our main sample compared between the TNG100DM (purple curve) and FIREboxDM
(black curves).