# Preliminary Evidence for Lensing-Induced Alignments of High-Redshift
Galaxies in JWST-CEERS

Viraj Pandya Hubble Fellow  Columbia Astrophysics Laboratory, Columbia
University, 550 West 120th Street, New York, NY 10027, USA  Abraham Loeb
Center for Astrophysics $|$ Harvard $\&$ Smithsonian, Cambridge, MA 02138, USA
Black Hole Initiative, Harvard University, Cambridge, MA 02138, USA  Elizabeth
J. McGrath Department of Physics and Astronomy, Colby College, Waterville, ME
04901, USA  Guillermo Barro Department of Physics, University of the Pacific,
Stockton, CA 90340 USA  Steven L. Finkelstein Department of Astronomy, The
University of Texas at Austin, Austin, TX, USA  Henry C. Ferguson Space
Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA
Norman A. Grogin Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218, USA  Jeyhan S. Kartaltepe Laboratory for Multiwavelength
Astrophysics, School of Physics and Astronomy, Rochester Institute of
Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA  Anton M.
Koekemoer Space Telescope Science Institute, 3700 San Martin Dr., Baltimore,
MD 21218, USA  Casey Papovich Department of Physics and Astronomy, Texas A&M
University, College Station, TX 77843-4242, USA  George P. and Cynthia Woods
Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M
University, College Station, TX 77843-4242, USA  Nor Pirzkal Space Telescope
Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  L. Y. Aaron
Yung Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD
21218, USA

###### Abstract

The majority of low-mass ($\log_{10}M_{*}/M_{\odot}=9-10$) galaxies at high
redshift ($z>1$) appear elongated in projection. We use JWST-CEERS
observations to explore the role of gravitational lensing in this puzzle. The
typical galaxy-galaxy lensing shear $\gamma\sim 1\%$ is too low to explain the
predominance of elongated early galaxies with ellipticity $e\approx 0.6$.
However, non-parametric quantile regression with Bayesian Additive Regression
Trees reveals hints of an excess of tangentially-aligned source-lens pairs
with $\gamma>10\%$. On larger scales, we also find evidence for weak lensing
shear. We rule out the null hypothesis of randomly oriented galaxies at
$\gtrsim 99\%$ significance in multiple NIRCam chips, modules and pointings.
The number of such regions is small and attributable to chance, but coherent
alignment patterns suggest otherwise. On the chip scale, the average complex
ellipticity $\langle e\rangle\sim 10\%$ is non-negligible and beyond the level
of our PSF uncertainties. The shear variance
$\langle\overline{\gamma}^{2}\rangle\sim 10^{-3}$ is an order of magnitude
above the conventional weak lensing regime but is more sensitive to PSF
systematics, intrinsic alignments, cosmic variance and other biases. Taking it
as an upper limit, the maximum implied “cosmic shear” is only a few percent
and cannot explain the elongated shapes of early galaxies. The alignments
themselves may arise from lensing by a protocluster or filament at $z\sim
0.75$ where we find an overabundance of massive lens galaxies. We recommend a
weak lensing search for overdensities in “blank” deep fields with JWST and the
Roman Space Telescope.

##  1 Introduction

Almost thirty years ago, Hubble Space Telescope (HST) observations revealed
that faint early galaxies preferentially appear elongated in projection (Cowie
et al., 1995; van den Bergh et al., 1996). This has since been confirmed with
newer instruments and statistical samples from larger surveys with HST
(Elmegreen et al., 2005; Ravindranath et al., 2006; van der Wel et al., 2014;
Zhang et al., 2019). Most recently, Pandya et al. (2024) revived interest in
this puzzle by showing that low-mass galaxies with
$\log_{10}M_{*}/M_{\odot}\sim 9-10$ at $z>1$ continue to appear preferentially
elongated with axis ratios $b/a\sim 0.3-0.6$ even in deeper observations with
the James Webb Space Telescope (JWST) at rest-frame optical wavelengths (see
also Kartaltepe et al., 2023; Robertson et al., 2023; Vega-Ferrero et al.,
2024).

These early elongated galaxies have continued to defy a clear explanation. The
safest bet is that there is a surface brightness selection effect against
detecting rounder face-on disks (Dalcanton & Shectman, 1996; Elmegreen et al.,
2005; Loeb, 2024) or that we are preferentially seeing an elongated star-
forming “proto-bar” but not the extended stellar disk. Zhang et al. (2019)
showed that at most $\sim 20\%$ of elongated galaxies would be missed by HST
if their light profiles were re-projected into the face-on view, which is not
enough to explain the high inferred elongated fractions of $\sim 50-70\%$.
Pandya et al. (2024) demonstrated that JWST should be even more complete to
disks with all orientations at $z\sim 1-8$ down to
$\log_{10}(M_{*}/M_{\odot})\sim 9$ over a reasonable range of sizes, but
stressed the need to scrutinize even deeper surveys. Alternatively, we may be
seeing tidal debris from ongoing interactions since merger rates are expected
to increase towards high redshift (Conselice, 2014; Somerville & Davé, 2015).
Finally, these objects may truly have prolate or triaxial 3D shapes that
reflect their formation via mergers along cosmic web filaments (Ceverino et
al., 2015; Tomassetti et al., 2016).

Here we explore a possibility not previously considered in the literature: are
early galaxies preferentially elongated because they have a higher probability
of being gravitationally lensed? It is already well known that the optical
depth to lensing increases for higher redshift systems (e.g., Barkana & Loeb,
2000) and that this results in a “magnification bias” which modifies the
number counts of high-redshift galaxies (Turner et al., 1984; Wyithe et al.,
2011; Mason et al., 2015). The larger covering fraction of possible foreground
lenses combined with bigger angular diameter distances to high-redshift
sources can also help maximize the Einstein radius for strong lensing (Treu,
2010). JWST has already identified multiple strongly lensed background
galaxies in massive clusters missed by Hubble (e.g., Pascale et al., 2022;
Mowla et al., 2022; Frye et al., 2023; Mowla et al., 2024; Bradley et al.,
2024), with roughly $\sim 25-65$ cases expected in “blank” JWST deep fields
(Holloway et al., 2023; Casey et al., 2023). One new Einstein ring has already
been detected in COSMOS-Web based on serendipitous visual inspection (Mercier
et al., 2023; van Dokkum et al., 2024).

Perhaps the most relevant connection is to weak lensing which is sensitive to
the assumed intrinsic shapes of galaxies. JWST is uniquely enabling the use of
lower mass, higher redshift sources for weak lensing studies around massive
clusters (Finner et al., 2023a; Harvey & Massey, 2024) but the intrinsic
shapes and intrinsic alignments of these preferentially elongated background
galaxies remain poorly understood. Pandya et al. (2019) proposed that if most
early galaxies are indeed nearly prolate and forming along cosmic web
filaments, then they are expected to show very strong intrinsic alignments
which may be an under-appreciated source of bias for weak lensing. On the
other hand, lensing by foreground large-scale structure may itself contribute
a non-negligible amount of “cosmic shear” that leads to some net elongation
and alignments of distant sources. Of course, the fact that it is the lower
mass systems at high-redshift that preferentially appear elongated suggests
that this phenomenon is due to their intrinsic shapes rather than lensing, but
it is still important to quantify the magnitude of the effect and possible
implications.

With these fundamental questions in mind, here we pursue a systematic search
for both galaxy-galaxy lensing candidates and large-scale alignments in
“blank” JWST deep fields. This is complementary to traditional lensing studies
around bright foreground clusters and involves looking for the statistical
correlations between ellipticity, orientation and shear that are the hallmarks
of gravitational lensing (e.g., Tyson et al., 1990; Brainerd et al., 1996).
For this pilot study, we will use the JWST Cosmic Evolution Early Release
Science survey (CEERS; Program ID 1345; Finkelstein et al., 2023), which was
selected to not have an obvious bright cluster in the foreground but may still
have overdensities at higher redshift. By averaging over the orientations of
background galaxies on multiple scales, we will demonstrate the potential of
JWST for constraining the presence of any such foreground mass concentration
even if it is “dark” (e.g., Kaiser et al., 2000; Bacon et al., 2000; Wittman
et al., 2000; Maoli et al., 2001; Rhodes et al., 2001; Refregier, 2003;
Kilbinger, 2015). With that said, it is not our goal to measure the precise
amplitude of any lensing signal but rather to place an upper limit on its
contribution to the elongation of early galaxies.

This paper is organized as follows. In section 2 we describe the data and in
section 3 we detail our methods. We present our results on galaxy-galaxy
lensing in section 4 and on large-scale alignments in section 5. We discuss
possible explanations and implications in 6. Finally, we summarize in section
7. We assume a standard Planck Collaboration et al. (2016) cosmology
throughout with $h=0.6774$, $\Omega_{\rm m,0}=0.3075$,
$\Omega_{\rm\Lambda,0}=0.691$ and $\Omega_{\rm b,0}=0.0486$.

##  2 Data

We use data from the JWST-CEERS survey (Program ID 1345; Finkelstein et al.,
2023).111These data can be found on MAST: http://dx.doi.org/10.17909/z7p0-8481
(catalog 10.17909/z7p0-8481) CEERS covers a portion of the Extended Groth
Strip (EGS; Groth et al., 1994; Davis et al., 2007) with NIRCam imaging in 10
pointings. Data reduction details are given in Bagley et al. (2023). Here we
use a source catalog derived by Finkelstein et al. (2023) with the original
SourceExtractor code (Bertin & Arnouts, 1996). We use photometric redshifts
and stellar masses derived from 13-filter spectral energy distribution (SED)
fitting by Barro et al. (2023) using the EAZY code (Brammer et al., 2008). The
13 filters include six broadband ones from NIRCam (F115W, F150W, F200W, F277W,
F356W, F444W), one medium-band NIRCam filter (F410M), six broadband filters
from HST-ACS/WFC3 (F606W, F814W, F105W, F125W, F140W, F160W). The multi-
wavelength SEDs give reasonably well constrained photometric redshifts and
stellar masses for our selected sample as described below in subsection 3.6.

Our galaxy shape measurements are based on single-component Sérsic fits with
galfit (Peng et al., 2002) by McGrath et al. (in prep.). This provides the key
quantities that we need for our lensing analysis: effective (half-light)
radius, projected axis ratio, and position angle of the major axis measured
for each galaxy independently in all six broadband NIRCam filters. The shape
catalog also includes empirical errors for these quantities that are derived
by matching each observed galaxy to 100 simulated sources which were inserted
into blank regions of the mosaic and which have similar magnitude, size and
Sérsic index. These empirical errors help quantify the systematic
uncertainties that generally dominate galaxy shape measurements and are not
accounted for by the formal statistical uncertainties from galfit. During the
Sérsic model fitting process, an empirical, filter-dependent global PSF was
used to recover “intrinsic” galaxy sizes, axis ratios and position angles
before convolution with a global empirical PSF. These global PSFs were created
by stacking stars throughout the CEERS footprint (see section 3.2 Finkelstein
et al., 2023). In Appendix A, we investigate the spatial dependence of the PSF
in CEERS by measuring the quadrupole moments of individual stars and quantify
the level of bias expected.

##  3 Methods

Here we describe our methods to study galaxy-galaxy lensing and large-scale
alignments (in that order).

###  3.1 Lens and Source Selection

We select possible foreground massive lens galaxies by requiring $z<1$,
$\log_{10}(M_{*}/M_{\odot})>10.5$ and good Sérsic fits (galfit flag of zero or
one) in the F115W filter which traces the rest-frame optical/near-infrared. Of
the 53885 galaxies available in the catalog, this leads to 77 massive lens
galaxies. Based on visual inspection, we discard one that is clearly a
deblending artifact resulting in 76 lenses.

For background source galaxies, we impose $z>1$ and
$\log_{10}(M_{*}/M_{\odot})>9$ which is above the CEERS galaxy completeness
limit to at least $z=8$ for a reasonable range of sizes and apparent
magnitudes (Appendix B of Pandya et al., 2024). We do not make any cut on
color or star formation rate. Of the 53385 sources in the catalog, this yields
6648 galaxies. We further require that the galaxies have good Sérsic fits
(galfit flag of zero or one) in the filter that most closely tracks the rest-
frame optical ($\sim 5000$Å) at their redshift. This means F115W for
$z=1.0-1.5$, F150W for $z=1.5-2.0$, F200W for $z=2-3$, F356W for $z=3-6$ and
F444W for $z>6$. We only include galaxies whose intrinsic (i.e., PSF-
deconvolved) Sérsic effective radius is greater than the PSF FWHM of their
assigned rest-optical filter. These cuts yield 4135 source galaxies but we
discard 267 that are clearly artifacts based on visual inspection, and another
20 that have catastrophically high uncertainty estimates of $\Delta(b/a)>1$,
$\Delta(\rm{PA})>30^{\circ}$ or $\Delta(r_{\rm e})>1^{\prime\prime}$. This
leaves 3848 sources of which 3073 ($\sim 80\%$) are low-mass with
$\log_{10}(M_{*}/M_{\odot})=9-10$.

###  3.2 Lens-Source Pair Selection

For every background galaxy, we start by computing its on-sky separation to
all 76 possible foreground lenses. We define a subset of these as “nearest”
source-lens pairs by assigning the closest on-sky lens to each background
galaxy. We verified that our results would be identical if we had instead
assigned to each source the lens with the highest predicted shear.

Figure 1 shows the demographics of our source-lens pairs in terms of the joint
distribution of lens redshift, source redshift and on-sky pair separation. The
lenses are roughly uniformly distributed over $z\sim 0.3-1.0$ but there
appears to be an overdensity at $z\sim 0.75$ (top gray histogram). The
background sources span a range of redshifts but drop off steeply from $z\sim
1-10$ (top cyan histogram). Interestingly, the highest redshift background
galaxies (orange/red points) are preferentially found near lenses with $z\sim
0.75$. It is unclear if this is due to magnification bias, a selection effect,
cosmic variance or some combination thereof. At a given lens redshift,
background sources span a range of on-sky pair separations $\theta\lesssim
250$ arcsec. For the average background galaxy, the nearest massive foreground
lens is $\sim 1$ arcmin away.

Figure 1: Joint distribution of lens redshift, source redshift and pair
separation for all background galaxies paired with their nearest on-sky
massive foreground lens. The lens redshift distribution is roughly uniform
between $z\sim 0.3-1.0$ with an overdensity at $z\sim 0.75$ (top gray
histogram). The source redshift distribution drops off rapidly at $z>1$ (top
cyan histogram). At fixed lens redshift, background galaxies span a range of
source redshifts and pair separations with $\theta\lesssim 250$ arcsec. On
average, background galaxies are $\sim 1$ arcmin from their nearest on-sky
massive foreground lens (right gray histogram).

###  3.3 Lens Model

We assume a singular isothermal sphere (SIS) lens model for all 76 massive
foreground galaxies. For this initial exploratory study, we believe this
strikes a balance between the simplicity of point lens models and the complex
shear analysis that is standard in weak lensing (Schneider et al., 1992).

Since we are interested in calculating the galaxy-galaxy lensing shear, we
begin by computing the angular diameter distances to the lens ($D_{\rm L}$),
to the source ($D_{\rm S}$), and between the lens and source ($D_{\rm LS}$).
Then we can compute the Einstein radius as (e.g., see Narayan & Bartelmann,
1996; Treu, 2010; Schneider, 2015):

$\theta_{\rm E}=4\pi\frac{D_{\rm LS}}{D_{\rm S}}\left(\frac{\sigma_{\rm
SIS}}{c}\right)^{2}$ (1)

where $c$ is the speed of light and $\sigma_{\rm SIS}$ is the (radially-
constant) 1D velocity dispersion of the SIS halo. For massive early-type
lenses, the central stellar velocity dispersion $\sigma_{*}$ has been shown to
correlate well with $\sigma_{\rm SIS}$ (i.e., being roughly equal; Bolton et
al., 2008; Zahid et al., 2016, 2018). We thus approximate the central stellar
velocity dispersion of all lenses assuming their centers are baryon-dominated
via

$\sigma_{*}=\sqrt{\frac{GM_{*}}{R_{\rm eff}^{\rm F115W}}}\;.$ (2)

This yields a roughly log-normal distribution of $\sigma_{*}$ with shape
parameter $s\sim 0.3$ and scale parameter $275$ km/s, which peaks at $\sim
250$ km/s and has a tail out to $\sim 580$ km/s. The distribution is
reasonable and compares well with, e.g., the range of $\sigma_{*}$ measured
for luminous red galaxies (LRGs) at $z<0.7$ in the SDSS-III/BOSS survey
(Thomas et al., 2013). In detail, since our lens sample extends down to
$\log_{10}M_{*}/M_{\odot}>10.5$, it likely includes many “fast rotator” early-
type galaxies (e.g., Cappellari, 2016). Equation 2 may need a correction to
account for the greater fraction of rotation in these stellar systems. We
forgo any such correction here, leaving it as a systematic that is only
relevant for our galaxy-galaxy lensing analysis but not cosmic shear
calculations (the latter is described next in subsection 3.5).

Taking $\sigma_{\rm SIS}=\sigma_{*}$, we can thus compute $\theta_{\rm E}$.
The shear experienced by a background galaxy a projected distance $\theta$
away from the center of a foreground lens is then simply

$\gamma=\frac{1}{2}\frac{\theta_{\rm E}}{\theta}\;.$ (3)

This shear can be thought of as a differential change in projected
ellipticity. It is larger for smaller on-sky source-lens separations and for
more massive lenses which have larger Einstein radii. For a given on-sky
separation $\theta$ and SIS lens, the shear is maximized as $D_{\rm LS}\to
D_{\rm S}$.

In addition to shear, we also compute the tangential alignment angle for every
source-lens pair as illustrated in Figure 2. Every background galaxy already
has a position angle $\phi$ for its major axis from the best-fitting Sérsic
model. We compute another position angle $\delta$ from the vector that
connects the centers of the background and lens galaxies. Both position angles
are defined with the same astronomical convention where $0^{\circ}$ is north
(up) and we limit the range between $-90^{\circ}$ and $90^{\circ}$ east of
north. The absolute difference between the two is $\psi\equiv|\phi-\delta|$
and lies between $0-90^{\circ}$. When $\psi\approx 0^{\circ}$, it means the
major axis of the background galaxy is pointing towards its nearest lens. When
$\psi\approx 90^{\circ}$, it instead means the major axis of the background
galaxy is perpendicular to the vector connecting the source-lens pair, i.e.,
we say the background galaxy is tangentially aligned with the lens.

Figure 2: Illustration of the various vectors and angles computed in this
work. We are primarily interested in the position angle $\phi$ of the major
axis of the background galaxy and its tangential alignment angle $\psi$ with
respect to the position of the lens. The on-sky angular pair separation is
$\theta$. All angles are defined between $-90^{\circ}$ to $90^{\circ}$ east of
north.

###  3.4 Bayesian Quantile Regression with BART

Since we are dealing with small sample sizes, we opt for a Bayesian analysis
of galaxy-galaxy lensing. Specifically, we will perform Bayesian “quantile
regression” which predicts the conditional distribution of different quantiles
of some property and their dependence on other variables. This provides more
information than linear regression which instead fits for a single mean
relation and is susceptible to outliers and non-linearities. Unlike simple
linear regression, quantile regression can also capture heteroscedasticity in
the data such as the scatter in $y$ changing with $x$. While quantile
regression can be done in a parametric way (e.g., fitting different lines for
different quantiles of $y$), here we opt for a non-parametric approach using
Bayesian Additive Regression Trees (BART; Chipman et al., 2010; Hill et al.,
2020; Martin et al., 2021). BART approximates a function by summing over many
small regression trees whose individual sizes and depths are subject to priors
that avoid overfitting.222In the limit of infinite regression trees, BART
approaches a Gaussian process but without the need to choose and tune a
covariance function (i.e., kernel). BART can also be thought of as a Bayesian
cousin of machine learning techniques like random forests and gradient-boosted
decision trees. In addition to the sum-of-trees, there is an error term that
is fit as part of the overall model.

We use the implementation of BART in the probabilistic programming language
PyMC (Quiroga et al., 2023; Oriol et al., 2023), which utilizes Hamiltonian
Monte Carlo for accelerated gradient-based sampling of the error term, and
particle Gibbs sampling for BART. BART has a number of hyper-parameters that
characterize its internal priors for the depths, splitting rules and values of
its trees to minimize overfitting. The default values have been shown to work
well for a variety of problems so we do not vary those. We use 100 trees to
fit our thousands of points in the $\log\gamma-e$ and $\log\gamma-\psi$
planes, but using 50 or 200 trees did not change our results. Note that
$\gamma$ has to be inferred while $e$ and $\psi$ are direct observables so we
treat the latter as “independent” variables.

We assume an asymmetric Laplace likelihood which is naturally well-suited for
quantile regression (Yu & Moyeed, 2001). It has three parameters: one
controlling the mean, another controlling the scale, and yet another
controlling the asymmetry. The location parameter of this likelihood is the
BART random variable with its default internal priors as described above. For
the scale parameter (which quantifies uncertainty around the BART-based mean),
we use a half-normal prior with standard deviation arbitrarily set to five.
Finally, the asymmetry parameter is set to the quantile we are trying to fit.
In this work, we will fit five quantiles: 0.1, 0.5, 0.9, 0.95 and 0.997.

We ran four chains with 3000 tuning (burn-in) and 3000 sampling draws. We
verified the convergence of the chains as follows. Trace plots revealed that
the posteriors of the scale parameter for the asymmetric Laplace likelihood
were similar from all four independent chains. Since the BART random variable
is really a collection of predictions from multiple trees for each of our
thousands of observed points, the recommended convergence check involves the
cumulative distribution function of the effective sample size and Gelman-Rubin
statistic for each chain. We verified that the effective sample size was $\gg
1000$ and the Gelman-Rubin statistic was $\lesssim 1.02$ for all BART
components of each chain, thus implying convergence.

###  3.5 Complex Ellipticity Analysis

We will average over the orientations of background galaxies both in bins of
galaxy-galaxy shear and more generally on the scale of individual NIRCam chips
($64^{\prime\prime}\times 64^{\prime\prime}$), modules ($2.2^{\prime}\times
2.2^{\prime}$) and pointings ($2.2^{\prime}\times 5.1^{\prime}$). For this, we
will follow standard practice in gravitational lensing studies (e.g.,
Schneider et al., 1992) and assign each galaxy a complex number

$e=|e|\exp(2i\phi)=|e|\cos(2\phi)+i|e|\sin(2\phi)$ (4)

where $|e|\equiv(1-b/a)$ is the usual projected ellipticity and $\phi$ is the
position angle measured east of north (up).333In weak lensing, it is common
practice to instead define $|e|=(a-b)/(a+b)$ or
$|e|=(a^{2}-b^{2})/(a^{2}+b^{2})$. Our results are not sensitive to the choice
of a particular definition for the norm $|e|$. The use of $2\phi$ accounts for
the symmetry of rotating an ellipse by $\pi$, i.e., galaxy shape vectors do
not point in one or the other direction for a given orientation. This
formalism makes it easy to compute the average ellipticity as the magnitude of
the average complex number $\langle e\rangle$ and the average orientation as
its phase.444It is common practice to divide individual ellipticities by a
“shear polarizability tensor” $P^{\gamma}$ or “shear responsivity factor”
$\mathcal{R}$ (Massey et al., 2007a). These are of order unity and may depend
on galaxy morphology so we neglect them here for simplicity.

We will also compute the “shear variance”
$\langle\overline{\gamma}^{2}\rangle$ which quantifies the degree of
correlation between the complex ellipticities of all possible pairs of
galaxies in some region. This summary statistic is frequently used in cosmic
shear studies (Refregier, 2003; Kilbinger, 2015) and can be thought of as a
simpler alternative to more sophisticated shear-shear correlation functions,
which otherwise depend continuously on pair separation. Specifically, we adapt
equation (8) of Hämmerle et al. (2002):

$\overline{\gamma}^{2}_{\rm reg}=\frac{1}{N_{\rm gal}(N_{\rm
gal}-1)}\sum_{i=1}^{N_{\rm gal}}\sum_{j=i+1}^{N_{\rm gal}}e_{i}e_{j}^{*}$ (5)

where the normalization accounts for the number of unique pairs. Note that the
inner summation is staggered by one to prevent double counting pairs and that
the product involves the complex conjugate $e_{j}^{*}$. This
$\overline{\gamma}^{2}_{\rm reg}$ is the shear variance for a single chip,
module or pointing, but we want the average over all regions of a given type:

$\langle\overline{\gamma}^{2}\rangle=\frac{1}{N_{\rm reg}}\sum_{n=1}^{N_{\rm
reg}}|\overline{\gamma}^{2}_{\rm reg,n}|$ (6)

where the sum involves the norm of $\overline{\gamma}^{2}_{\rm reg,n}$, which
in general is a complex number. We will compute this shear variance for both
galaxies and our PSF stars (Appendix A) with the latter serving as null tests
for systematics.

###  3.6 Systematics Checks and Error Propagation

We can adapt the shear variance calculation of equation 5 to compute the star-
galaxy cross-correlation as a check of systematics:

$e_{\rm gal}\star e_{\rm star}^{*}=\frac{1}{N_{\rm gal}N_{\rm
star}}\sum_{i=1}^{N_{\rm gal}}\sum_{j=1}^{N_{\rm
star}}e_{\rm{gal},i}e_{\rm{star},j}^{*}\;.$ (7)

For this, we measured the complex ellipticity of individual stars throughout
CEERS using quadrupole moments as described in Appendix A. This provides an
important null test for baselining possible PSF systematics but we caution
that we only have $\sim 130$ stars throughout the survey footprint and
typically only $\sim 1-2$ stars per NIRCam chip, with $\sim 30\%$ of chips
having no star. To ensure that galaxies are only paired up with stars in their
corresponding rest-frame optical filter, we will compute this separately for
each NIRCam filter. Since we want to know the fractional contribution of the
star-galaxy cross-correlation to the total shear variance, we record the ratio
of the norms $|e_{\rm gal}\star e_{\rm star}^{*}|/|\overline{\gamma}^{2}_{\rm
reg}|$ for all individual regions before taking the average as in 6. This
approach of averaging ratios should be more robust especially on smaller
scales where “shot noise” from the galaxy ellipticity distribution can
dominate the variance between regions of the same size. Note that our galaxy
shapes were measured by fitting a Sérsic model convolved with a global
empirical PSF and that we do not attempt to perform local PSF corrections here
(but see Appendix A). The star-galaxy cross-correlation thus reflects any
additional local PSF contamination unaccounted for by the global PSF.

We use an efficient vectorized Monte Carlo approach to propagate errors on
$\langle e\rangle$, $\langle\overline{\gamma}^{2}\rangle$ and the star-galaxy
cross-correlation. When computing these summary statistics for any sample of
$N$ galaxies or stars, we first create 1000 random realizations of size $N$
each. For every observed galaxy, we randomly draw an ellipticity, position
angle and any other quantity of interest from a Gaussian with mean equal to
the fiducial catalog value and standard deviation equal to the empirical error
from McGrath et al. (in prep.). For stars, we do the same thing but in lieu of
errors on quadrupole moments, we use the standard deviation of $e_{1}$ and
$e_{2}$ from the star sample itself as a measure of uncertainty on those
quantities. For simplicity, we do not propagate errors on galaxy redshift or
mass since our analysis is not done in fine $z$ or $M_{*}$ intervals.555The
typical redshift uncertainty for our sources is $\Delta(z_{\rm
phot})/(1+z_{\rm phot})\approx 0.046\pm 0.082$ where $\Delta(z)$ is the
difference between the 84th and 16th percentile of the photometric redshift
posterior for each galaxy. Then we can compute our summary statistics in all
realizations and take the mean and standard deviation with the latter
quantifying the uncertainty due to galaxy/star shape error propagation. It is
almost always the case that the summary statistics are well constrained, e.g.,
$\mu_{\langle e\rangle}/\sigma_{\langle e\rangle}\gg 1$. However, this cannot,
by itself, be used to assess the significance of an alignment signal. For
that, we need a null hypothesis test to rule out randomly oriented galaxies,
which we describe in subsection 3.7.

Our Monte Carlo error propagation method gives a formal uncertainty for
summary statistics computed on the scale of the entire survey. However, we
only have one survey limited to a small ($\sim 0.028$ deg2) part of the sky so
our shear variance $\langle\overline{\gamma}^{2}\rangle$ is almost certainly
dominated by cosmic variance. It is not our goal in this paper to constrain
cosmological parameters or perform detailed comparisons to theory so we do not
attempt to estimate the cosmic variance uncertainty on
$\langle\overline{\gamma}^{2}\rangle$. Instead, we will simply place an upper
limit on the weak lensing shear experienced by high-redshift galaxies. We
defer an estimation of the cosmic variance uncertainty to future weak lensing
analyses which can empirically constrain the field-to-field variance using
multiple “blank” JWST deep fields.

###  3.7 Null Hypothesis Test for Alignments

For any sample of galaxies that show elevated $\langle e\rangle$ suggestive of
alignments, we need to compare their observed $\langle e\rangle$ to the
distribution of $\langle e\rangle$ expected under the null hypothesis that
they are randomly oriented. For $N$ observed galaxies, we create 1 million
realizations, each of which is assigned $N$ random complex ellipticity
vectors. We randomly draw orientations $\phi$ from a uniform distribution
$\mathcal{U}(-90^{\circ},90^{\circ})$ which should “break” any alignments (if
present). Since there is no good, universal distribution from which to
randomly draw ellipticities, we fix the distribution of complex magnitudes $e$
to the actual observed ellipticities of the galaxies.666A uniform distribution
is clearly not a good choice for high-redshift, low-mass bins which are biased
towards high $e$ (Pandya et al., 2024) or for high-mass spheroid-dominated
bins which are biased towards low $e$ (e.g. Chang et al., 2013). This lets us
compute the distribution of the average $\langle e\rangle$ under the null
hypothesis that the $N$ galaxies are randomly oriented. We then calculate the
$p$-value for a given chip as the fraction of its random realizations where
the null $\langle e\rangle$ exceeds the observed $\langle e\rangle$. This
procedure is used both for our galaxy-galaxy lensing candidates and when
averaging over the orientations of all background galaxies in a given NIRCam
chip, module, pointing and the scale of the entire survey.

##  4 Galaxy-Galaxy Lensing

In this section we present our results on galaxy-galaxy lensing in JWST-CEERS.

###  4.1 Shear versus ellipticity

Figure 3 shows the joint distribution of galaxy-galaxy lensing shear777Recall
that $\gamma=0.5\theta_{\rm E}/\theta$. The distribution of $\theta_{\rm E}$
itself for all source-lens pairs is roughly a log-normal with mean $\mu\sim
0.4^{\prime\prime}$ and scale $\sigma\sim 1.1$. It peaks at $\theta_{\rm
E}\approx 0.7^{\prime\prime}$ and extends to $7^{\prime\prime}$. and source
ellipticity. The ellipticity distribution of background galaxies is clearly
non-uniform and biased towards elongated objects with $e\sim 0.6$ on average.
The ellipticity distribution of massive sources alone is flatter, consistent
with them being randomly oriented disks, but this implies that the low-mass
sources are the preferentially elongated ones. Most pairs have very small
shear with the mean being $\langle\gamma\rangle\approx 0.02$ consistent with
what is expected for weak lensing (e.g., Oguri et al., 2012). Because the mean
ellipticity of background galaxies is much larger than the typical expected
shear, galaxy-galaxy lensing cannot be the primary driver for the preferential
elongation of low-mass high-redshift galaxies.

However, there is a non-negligible tail towards larger shears in excess of
$0.1$ which is beyond the conventional weak lensing regime. One object that is
highly elongated with $e\approx 0.7$ has an unusually large shear of
$\gamma\approx 0.6$. Our Bayesian non-parametric quantile regression with BART
reveals that shear does not correlate with ellipticity for the majority of
pairs below the 95% quantile of shear. Even for shears in the $99.7$%
quantile, the data does not show a strong correlation with ellipticity and the
Bayesian credible interval is large. At very high $e>0.7$, we see a downturn
in the median $\gamma$ of the 99.7% shear quantile, but the uncertainty is
large due to the small number of such extreme objects. We confirmed that the
$19$ objects with $e>0.85$ and $\gamma<0.1$ are in fact highly elongated
galaxies and not artifacts (coincidentally, two of them are shown as prolate
candidates in Figure 17 of Pandya et al., 2024). The tail of large shear pairs
motivates looking for correlations with the source-lens tangential alignment
angle, which we turn to next.

Figure 3: Joint distribution of shear and ellipticity for all nearest source-
lens pairs (black points). The colored solid lines show the conditional
distribution of shear on ellipticity in five quantiles from our non-parametric
quantile regression with Bayesian Additive Regression Trees: 10% (blue), 50%
(orange), 90% (green), 95% (magenta) and 99.7% (red). The light and dark
shaded regions reflect the 50% and 95% credible intervals, respectively. Most
background galaxies are predicted to have very small shear but there is a
strong tail towards $\gamma>0.1$. The data does not show a strong correlation
between shear and ellipticity regardless of shear quantile. The projected
ellipticity distribution of background galaxies is non-uniform and biased
towards high values (top gray histogram). If we restrict to high-mass
background galaxies, they show a flatter distribution (top cyan histogram) as
expected for randomly oriented disks, meaning that it is the low-mass galaxies
that are preferentially elongated.

###  4.2 Shear versus tangential alignment angle

Figure 4 plots the joint distribution of shear $\gamma$ and tangential
alignment angle $\psi$. Most pairs follow a broad, relatively uniform
distribution of $\psi$ (top gray histogram). There is a hint that high-shear
pairs with $\gamma>0.2$ tend to be tangentially aligned with larger $\psi$,
but this is only significant at the $\sim 96.6\%$ level ($p=0.0337$) based on
a two-sample Kolmogorov-Smirnov test using the unbinned, unweighted
distributions of $\psi$. While this hypothesis test based on the 1D
marginalized distribution of $\psi$ is useful, it is not entirely appropriate
since we are dealing with multi-dimensional data and really looking for
correlations (Feigelson & Babu, 2012).

Our Bayesian non-parametric quantile regression confirms that most pairs have
low shear (the 95% quantile is $\gamma\lesssim 0.1$), and that those low-shear
pairs have no correlation with alignment angle. However, the conditional
distribution of the extreme $99.7$% quantile of shear shows hints of a
positive correlation with alignment angle although our Bayesian credible
intervals (i.e., uncertainties) are large.888This is mainly driven by one
“outlier” with $\gamma\sim 0.38$ and $\psi\sim 0.5^{\circ}$. Visual inspection
reveals that this is clearly a high-redshift dropout galaxy that only appears
in redder NIRCam filters and is almost perfectly radially-aligned with a
massive foreground galaxy $\sim 2.8^{\prime\prime}$ away. For pairs in the
$99.7\%$ quantile of shear, when $\psi$ is large, the shear tends to be larger
by $\sim 5\%$. There is also a hint of heteroscedasticity in the data such
that the scatter in shear increases with $\psi$. In other words, at small
$\psi$, the 99.7% quantile of shear may simply be a tail of the log-normal
distribution of shear of all pairs, but at large $\psi$, the extremes of the
conditional shear distribution may comprise a broader tail. This would imply
that among the many tangentially-aligned pairs, there may be some that are
lensed, precisely those with high shear. With that said, our uncertainties are
large and this may be due to both our small sample sizes and lensing
systematics. By combining data from multiple “blank” JWST deep fields, future
analyses can adapt our approach to help tighten constraints on a possible
excess of tangentially-aligned source-lens pairs.

Figure 4: Joint distribution of shear and tangential alignment angle for all
nearest source-lens pairs (black points). The colored lines show the
conditional distribution of shear on tangential alignment angle in five
quantiles from our non-parametric quantile regression with Bayesian Additive
Regression Trees: 10% (blue), 50% (orange), 90% (green), 95% (magenta) and
99.7% (red). The light and dark shaded regions reflect the 50% and 95%
Bayesian credible intervals, respectively. Most pairs have small shear (right
gray histogram) and there is no correlation with $\psi$ even in the $95\%$
quantile of shear (magenta curve). In contrast, there is a hint of an excess
of tangentially-aligned pairs at high shear with $\gamma>0.2$ (top cyan
histogram). There is also a hint of a positive correlation in the 99.7% shear
quantile (red curve) wherein $\gamma$ tends to be $\sim 5\%$ higher for
tangentially-aligned pairs but the uncertainties are large. We will refer to
pairs with $\gamma>0.15$ and $\psi>75^{\circ}$ as tangentially-aligned lensing
candidates.

###  4.3 Images of lensing candidates

Figure 5 shows images of the 17 source-lens pairs that have $\gamma>0.1$ and
$\psi>75^{\circ}$. We will refer to these as our tangentially-aligned galaxy-
galaxy lensing candidates. Table 1 lists their properties to facilitate
follow-up observations. The cut on $\gamma>0.1$ places these candidates in a
shear regime that is an order of magnitude above the conventional weak lensing
limit. Many of them are elongated and some even show arc-like distortions that
may be suggestive of intermediate lensing. None of these are truly in the
strong lensing regime (i.e., all have $\theta>\theta_{\rm E}$) except possibly
one, which we now discuss.

The highest shear candidate is a $\sim 10^{9.5}M_{\odot}$ galaxy that is
clearly elongated with $e\sim 0.72$ and has a photometric redshift $z\sim
2.1$. It is nearly perfectly tangentially aligned ($\psi\approx 88^{\circ}$)
and lies within the Einstein radius of its associated lens ($\theta\sim
3.75^{\prime\prime}$, $\theta_{\rm E}\sim 4.75^{\prime\prime}$), making it our
only strong lensing candidate within the JWST-CEERS footprint. If its
associated lens, which happens to be our most massive one with the highest
$\sigma_{*}$, satisfies our SIS assumption, then we expect another brighter
image on the opposite side outside the Einstein radius with an image
separation $\Delta\theta=2\theta_{\rm E}=9.5^{\prime\prime}$. We do not find
any obvious source there with a similar redshift as the main image. The
massive foreground lens galaxy itself is spectroscopically-confirmed to be at
$z\sim 0.78$. If we overestimated its stellar velocity dispersion by a factor
of two, then with all else fixed, $\theta_{\rm E}$ would be halved and the
source would no longer lie within the Einstein radius so a second image would
not be expected (it would still have $\gamma\sim 0.3$ in this case). Of
course, another possibility is that the photometric redshift of our lensing
candidate is very wrong and it is instead physically associated with the lens
as a satellite. Note that the rest of our results in this paper do not hinge
on whether or not this candidate is confirmed (namely, the large-scale
alignments in section 5 are independent of this galaxy-galaxy lensing issue).

Figure 5: Example $3^{\prime\prime}\times 3^{\prime\prime}$ false-color RGB (F115W+F200W+F356W) postage stamps of our 17 tangentially-aligned lensing candidates with $\gamma>0.1$ and $\psi>75^{\circ}$. The gray arrows point in the direction of their closest foreground massive lens (typically beyond the cutout) and cyan arrows denote the major axis position angle. The large panel shows a $10^{\prime\prime}\times 10^{\prime\prime}$ cutout around our one strong lensing candidate. The solid white circle shows our fiducial Einstein radius with the source clearly inside it but there is no obvious second image on the opposite (lower-right) side of the lens. If the Einstein radius was halved (dashed white circle), a second image would not be expected and the shear would drop to $\gamma\sim 0.3$. Source ID | Source RA | Source Dec |  Source $z$ |  Source $\log_{10}(M_{*}/M_{\odot})$ |  Source $R_{\rm e}$ |  Source $e$ |  Source $\phi$ | Rest Filter | $\theta$ |  |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
23392 | 214.80777 | 52.86992 | 2.619 | 9.07 | 1.39 | 0.67 | -69.17 | F200W | 3.64 |  |  |  |  |  |  |  |  |  |
24257 | 214.8167 | 52.86969 | 1.118 | 9.43 | 1.9 | 0.14 | 78.41 | F115W | 2.61 |  |  |  |  |  |  |  |  |  |
24459 | 214.81947 | 52.87288 | 2.338 | 11.23 | 6.67 | 0.4 | -43.11 | F200W | 13.73 |  |  |  |  |  |  |  |  |  |
24638 | 214.83818 | 52.88739 | 2.359 | 10.89 | 1.14 | 0.38 | -64.52 | F200W | 16.12 |  |  |  |  |  |  |  |  |  |
28283 | 214.85403 | 52.87866 | 3.382 | 9.45 | 1.6 | 0.81 | -48.63 | F356W | 10.4 |  |  |  |  |  |  |  |  |  |
29037 | 214.85188 | 52.87375 | 2.349 | 9.06 | 1.04 | 0.58 | 81.11 | F200W | 11.07 |  |  |  |  |  |  |  |  |  |
51412 | 214.98434 | 52.89845 | 1.805 | 10.23 | 3.3 | 0.49 | 34.95 | F150W | 4.74 |  |  |  |  |  |  |  |  |  |
54273 | 214.81218 | 52.81713 | 5.119 | 9.43 | 1.31 | 0.52 | 43.65 | F356W | 14.68 |  |  |  |  |  |  |  |  |  |
65925 | 215.1386 | 52.98904 | 2.311 | 9.22 | 0.7 | 0.51 | 61.17 | F200W | 5.69 |  |  |  |  |  |  |  |  |  |
70412 | 215.09525 | 52.93305 | 1.995 | 10.06 | 2.21 | 0.61 | -67.7 | F150W | 2.99 |  |  |  |  |  |  |  |  |  |
86621 | 214.86714 | 52.79953 | 2.121 | 9.47 | 5.9 | 0.72 | 60.04 | F200W | 3.75 |  |  |  |  |  |  |  |  |  |
86695 | 214.87056 | 52.80254 | 4.764 | 9.49 | 1.61 | 0.25 | -53.02 | F356W | 16.82 |  |  |  |  |  |  |  |  |  |
87402 | 214.87488 | 52.80144 | 2.85 | 10.27 | 2.43 | 0.37 | -28.9 | F200W | 27.16 |  |  |  |  |  |  |  |  |  |
94812 | 214.84508 | 52.79176 | 4.247 | 9.7 | 1.6 | 0.61 | 29.23 | F356W | 10.17 |  |  |  |  |  |  |  |  |  |
95439 | 214.86063 | 52.80203 | 1.654 | 9.07 | 1.91 | 0.75 | 7.4 | F150W | 5.75 |  |  |  |  |  |  |  |  |  |
95498 | 214.84607 | 52.78814 | 2.622 | 9.59 | 1.11 | 0.63 | -73.42 | F200W | 11.79 |  |  |  |  |  |  |  |  |  |
95846 | 214.86427 | 52.79873 | 1.316 | 10.66 | 3.49 | 0.36 | -4.83 | F115W | 12.8 |  |  |  |  |  |  |  |  |  |
Table 1: Properties of our 17 tangentially-aligned lensing candidates and
their main associated lenses. All IDs and positions are from the CEERS
photometric catalog with the latter in decimal degrees. Photometric redshifts
and stellar masses are from SED fitting, effective Sérsic radii are in proper
kpc, and the lens stellar velocity dispersions are in km/s. Position angles
are defined between $-90^{\circ}<\phi<90^{\circ}$ east of north, and alignment
angles are between $0^{\circ}<\psi<90^{\circ}$ with $0^{\circ}$ for radial and
$90^{\circ}$ for tangential alignment. Pair separations and Einstein radii are
in arcsec. The full table is available for download from the journal.

###  4.4 Source magnifications and demographics

Figure 6 shows that all of our lensing candidates fall within twice the
projected virial radius of their associated lens. Here we have estimated the
projected halo virial radius assuming $R_{\rm eff}/R_{\rm vir}=0.02$
(Kravtsov, 2013; Somerville et al., 2018). Importantly, the shear of our
lensing candidates falls off with “impact parameter” $\theta/R_{\rm vir}^{\rm
proj}$ as $\gamma\propto\theta^{-1}$, as it should. The scatter around this
relation must come from additional $\theta_{\rm E}$ dependence in Equation 3
or due to our simple scatter-free approximation that $R_{\rm eff}/R_{\rm
vir}=0.02$. We can estimate the magnification due to an SIS lens as

$\mu_{\pm}=\theta_{\pm}/\beta$ (8)

where $\beta=\theta_{\pm}-\theta_{\rm E}$ is the “true” source position,
$\theta_{\pm}$ are the two image positions and $\mu_{\pm}$ are their
corresponding magnifications (Narayan & Bartelmann, 1996). All of our sources,
except possibly the one strong lensing candidate, lie outside their Einstein
radii so there is only one image $\theta_{+}$ and these all have $\mu\sim
1-2$. Our strong lensing candidate has $\mu_{+}\sim-3.7$ for the inner image
(where the negative sign means the image parity is flipped with respect to the
true source) and the second image is predicted to be $\sim 1.5\times$ brighter
with $\mu_{-}\sim 5.7$.

Figure 6: Radial dependence of shear on “impact parameter” for all source-lens
pairs (small circles) and for our tangentially-aligned lensing candidates with
$\gamma>0.1$ and $\psi>75^{\circ}$ (stars). All of our lensing candidates fall
within twice the projected virial radius of their associated foreground lens
and their shear drops off as $\gamma\propto\theta^{-1}$. The colorbar denotes
the lensing magnification $\mu$ which is $\sim 1-2$ for most candidates. If
our highest shear object is truly strongly lensed, the inner image has
$|\mu|\sim 3.7$ and the outer image is predicted to have $\mu\sim 5.7$ (small
inset square).

Figure 7 shows the joint distribution of lens $\sigma_{*}$, the angular
diameter distance ratio $D_{\rm LS}/D_{\rm S}$ and the pair separation
$\theta$, all of which go into predicting the shear (equation 3). The lensing
candidates are associated with preferentially more massive lenses, smaller on-
sky separations and higher $D_{\rm LS}/D_{\rm S}$ since that maximizes the
Einstein radius $\theta_{\rm E}$ for an SIS lens with a given $\sigma_{*}$. In
other words, the distribution of $\theta_{\rm E}$ is biased towards larger
values ($\sim 1.5-6^{\prime\prime}$) for lensing candidates compared to the
overall log-normal distribution given in Footnote 7.

Figure 7 also plots the joint distribution of source redshift, source apparent
magnitude and source Sérsic effective radius with the latter two in the NIRCam
filter that most closely tracks the rest-frame optical. The lensing candidates
appear well mixed with the general population of sources except that they are
biased towards $z\gtrsim 2$ which just reflects the need for some minimal
source-lens angular diameter distance so that $\theta_{\rm E}$ is maximized.
The highest redshift candidate is at $z\sim 5$, and all candidates have
apparent magnitudes in the range $\sim 22-26$ AB mag in their rest-optical
filters. The candidates have rest-optical $R_{\rm eff}\sim 0.1-0.5$ arcsec
which requires adaptive optics or space-based telescopes for observational
follow-up.

Figure 7: Demographics of our lensing candidates relative to all source-lens
pairs. Left: Joint distribution of SIS lens velocity dispersion, angular
diameter distance ratio $D_{\rm LS}/D_{\rm S}$ and on-sky pair separation
$\theta$ for all source-lens pairs (small gray squares). These three
quantities set the level of shear. Darker colors correspond to smaller pair
separations. The large crosses are for our lensing candidates with
$\gamma>0.1$ and $\psi>60^{\circ}$. Lensing candidates are preferentially
found around more massive lenses (higher $\sigma_{*}$), at smaller pair
separations (darker colors), and at greater angular diameter distance from
their nearest lens which maximizes the Einstein radius $\theta_{\rm E}$.
Right: Joint distribution of source redshift, source apparent magnitude and
source Sérsic effective radius with the latter two measured in the NIRCam
filter that most closely tracks the rest-frame optical. Lensing candidates
roughly track the underlying population of source-lens pairs except they are
biased towards $z\gtrsim 2$ to maximize the angular diameter distance to our
lenses at $z<1$. Their rest-optical sizes are $\sim 0.1-0.5$ arcsec and
apparent magnitudes $\sim 22-26$ AB mag.

###  4.5 Correlated orientations of lensing candidates

The possible excess of tangentially-aligned source-lens pairs in the 99.7%
quantile of shear in Figure 4 motivates averaging over the orientations of
background galaxies in bins of galaxy-galaxy lensing shear. For galaxies with
negligible ($\gamma<0.01$) and weak ($0.01<\gamma<0.1$) shear, we expect the
magnitude of their average complex ellipticity vector to be negligible,
assuming they are randomly oriented and there are no systematics. Figure 8
confirms this even when splitting by mass. However, at moderate
($0.1<\gamma<0.2$) and high ($\gamma>0.2$) shear, we start to see evidence for
a non-zero average complex ellipticity and hence preferred orientation. The
standard errors come from our Monte Carlo error propagation method (subsection
3.6) and show that the trend is well constrained modulo any systematics. We
compare these observed averages to the null distribution of $\langle e\rangle$
assuming the galaxies are randomly oriented (subsection 3.7). Given the small
sample sizes, only the combined and low-mass subsamples in the moderate shear
bin remain statistically significant at the $>95\%$ level (the $p$-values are
reported above each point in Figure 8).

This follows the expectation that as shear increases, there must be some
residual orientation dictated by lensing. The complex average $\langle
e\rangle\sim 10-30\%$ is at least an order of magnitude above the conventional
weak lensing regime and unlikely to be due to PSF systematics (Appendix A).
These coherent alignments are surprising because our sources trace a number of
different foreground lens galaxies. Of course, if this were repeated using
galaxy-galaxy lensing candidates from multiple widely-spaced surveys, we would
expect $\langle e\rangle\to 0$ even in the moderate and strong shear regimes.
But the fact that our galaxy-galaxy lensing candidates have correlated
orientations within the relatively small ($\sim 0.028$ deg2) CEERS footprint
motivates a more detailed spatial analysis of alignments on larger scales,
which we turn to next.

Figure 8: Average ellipticity and orientation as a function of galaxy-galaxy
lensing shear. Left: Magnitude of the average complex ellipticity vector in
increasing bins of shear from left to right: negligible ($\gamma<0.01$), weak
($0.01<\gamma<0.1$), moderate ($0.1<\gamma<0.2$) and strong ($\gamma>0.2$).
All points reflect the mean and standard error from our Monte Carlo error
propagation method (subsection 3.6). Our combined background sample is shown
in black, high-mass only in red, low-mass only in blue, and $z>2$ low-mass in
green. At small shear, the average ellipticity is very small as it should be
for randomly oriented galaxies. At larger shear, we find a significantly
larger average residual ellipticity. The $p$-value from our null hypothesis
test is written above each point. Given the small sample sizes, only the two
moderate shear points enclosed in circles are significant at the $>95\%$
level. Right: Complex ellipticity vectors in bins of shear for individual
galaxies and the complex average of the combined sample (black lines). For
visualization, the symmetric complex ellipticity vectors have been “folded”
back onto the range $-90^{\circ}$ to $90^{\circ}$ measured east of north (up).
The magnitude of the average complex ellipticity vector is very small for the
negligible and weak shear bins, but is clearly non-zero with a preferred
orientation for the two higher shear bins.

##  5 Large-Scale Galaxy Alignments

In this section, we present results on large-scale alignments by averaging
over the complex ellipticities of background galaxies on multiple scales
defined by the NIRCam instrument. We start with individual NIRCam chips
($64^{\prime\prime}\times 64^{\prime\prime}$), then individual modules
($2.2^{\prime}\times 2.2^{\prime}$), followed by entire pointings
($2.2^{\prime}\times 5.1^{\prime}$ with a $\sim 0.7^{\prime}$ gap), and
finally the scale of the entire survey ($\sim 30^{\prime}\times 6^{\prime}$).
For this pilot study, restricting ourselves to the on-sky survey geometry
imposed by NIRCam mitigates systematics from detector gaps and otherwise
arbitrary binning. Appendix A quantifies the expected bias from PSF
systematics.

###  5.1 Alignments on NIRCam chip scales $(64^{\prime\prime}\times
64^{\prime\prime})$

Figure 9 shows that, on average, there are $\sim 50-100$ background galaxies
per arcmin2 that satisfy our $\log_{10}M_{*}/M_{\odot}>9$ cut in CEERS, which
covers a total area of about $\sim 100$ arcmin2 ($\sim 0.028$ deg2). This is
sufficient for weak lensing and cosmic shear studies. First, we show the on-
sky positions of our galaxy-galaxy lensing candidates with $\gamma>0.1$ and
$\psi>75^{\circ}$ compared to the foreground lenses and our overall background
source sample. The galaxy-galaxy lensing candidates appear clustered near the
southwest, which may help explain why we see correlated orientations in Figure
8. The lenses themselves appear roughly uniformly distributed throughout the
field, though we showed in Figure 1 that there is an overabundance of lenses
at $z\sim 0.75$. There is generally at least one, if not multiple, $z\sim
0.75$ foreground lenses near each lensing candidate. The projected $R_{\rm
vir}$ of our lenses are $\sim 0.1^{\prime}-3.5^{\prime}$, i.e., extending over
a single NIRCam chip or entire module.

Figure 9 also shows the magnitude of the average complex ellipticity when
averaging over background galaxies in individual NIRCam chips. The average
complex ellipticity is non-negligible and clearly in excess of $\langle
e\rangle\sim 10\%$ in many bins. This is an order of magnitude above the
conventional weak lensing regime and implies strong alignments in the
ellipticities and orientations of background galaxies. The standard error on
the mean $\langle e\rangle$ from our Monte Carlo error propagation is quite
low so in the next subsection we will use statistical null tests to quantify
the significance of these alignments.

Figure 10 splits the chip-scale “shear map” into three separate ones for high-
mass, low-mass, and high-redshift ($z>2$) low-mass background galaxies. There
are large number of NIRCam chips where high-mass galaxies show strong
alignments with $\langle e\rangle>10\%$. Despite the relative scarcity of
high-mass galaxies, the standard error on $\langle e\rangle$ from our Monte
Carlo error propagation is quite low, though we are likely dominated by
systematics. The low-mass sample down to $z>1$ shows non-negligible $\langle
e\rangle\gg 0.01$ similar to the combined map from Figure 9. Interestingly,
restricting to just low-mass galaxies at $z>2$ (in the predominantly elongated
regime; van der Wel et al., 2014; Zhang et al., 2019; Pandya et al., 2024)
reveals preferentially larger $\langle e\rangle\sim 20-30\%$ in multiple
regions compared to the $z<2$ low-mass sample. Whether this is telling us
something about lensing by a large foreground mass concentration, intrinsic
alignments or systematics requires us to first assess its statistical
significance, which we turn to next.

Figure 9: Projected on-sky lensing maps for CEERS in the standard “north-up”
frame. Left: Positions of the galaxy-galaxy lensing candidates (cyan pluses),
foreground lenses (circles colored by redshift) and our overall background
galaxy sample (grey histogram). The gray boxes show that there are $\sim
50-100$ per arcmin2 in each NIRCam chip ($64^{\prime\prime}\times
64^{\prime\prime}\sim 0.0003$ deg2) across this $\sim 100$ arcmin2 survey. The
lensing candidates, including the one possible strong lens (yellow star),
appear clustered near the southwest. The foreground lenses are roughly
uniformly distributed but Figure 1 suggests an overdensity at $z\sim 0.75$.
There is generally at least one $z\sim 0.75$ lens near each lensing candidate.
Right: Magnitude of the average complex ellipticity $\langle e\rangle$ when
averaging over the orientations of background galaxies in individual NIRCam
chips. The average complex ellipticity is non-negligible and reaches values in
excess of $10\%$ in some parts of the map. The standard error on the mean
(inset panel) is generally quite low but we are likely dominated by
systematics.

Figure 10: Same chip-scale “shear map” as in Figure 9 but now split into high-
mass (top), low-mass (middle), and high-redshift low-mass (bottom) background
galaxies. Only bins with $\geq 5$ galaxies are used. The top-right inset
panels show the number of galaxies arcmin-2 in each NIRCam chip. Top: Massive
sources show an elevated $\langle e\rangle\sim 0.3$ in many chips. Middle:
Low-mass sources at any $z>1$ also show non-negligible $\langle e\rangle\gg
0.01$ in many regions. Bottom: Low-mass sources at $z>2$ show larger $\langle
e\rangle\gtrsim 10\%$ in more chips compared to the $z<2$ low-mass sources.

###  5.2 Statistical significance from null hypothesis tests

Figure 11 quantifies the statistical significance of our results with the null
hypothesis test described in subsection 3.7. We find that the majority of
chips do not have statistically significant alignments since $p\gg 0.05$.
However, there are multiple cells where the null hypothesis can be ruled out
at $>95\%$ significance, including a handful at $>99\%$ significance. All of
these have $\langle e\rangle\gtrsim 0.1$ which, again, is at least an order of
magnitude larger than the conventional weak lensing regime. Interestingly,
some of these significant regions are in the southwest where we saw the
clustering of tangentially-aligned galaxy-galaxy lensing candidates in Figure
9.

In addition to “internal” chip-scale significance, it is useful to think about
the “map-scale” significance. Given that we have $80$ NIRCam chips, we expect
$\sim 80\times 0.05\approx 4$ detections by random chance at the $95\%$
confidence level assuming every bin is independent. We have more detections
than this simple threshold in the combined (6), high-mass (6) and $z>2$ low-
mass (8) maps. However, the low-mass map has only 3 detected chips and even
though all of these are at $\gtrsim 99\%$ significance, on the scale of the
map this may be a statistical fluke. Formally, when we use the Benjamini-
Hochberg false discovery rate and Bonferroni methods to “correct” the
individual chip-scale $p$-values to account for the large number (80) of tests
being done, none of the chips in any of the subsamples remain significant at
the $>95\%$ level. However, these conservative statistical methods may not be
using all of the available information (e.g., the clustering of chips with
detections) so it is also useful to look at images of individual chips with
alignments, which we turn to next.

Figure 11: Statistical significance of chip-scale “shear maps” for all
background galaxies combined (left column), high-mass sources (second-from-
left column), low-mass sources (second-from-right column), and $z>2$ low-mass
sources (right column). The top row shows the $p$-value in individual
$64^{\prime\prime}\times 64^{\prime\prime}$ chips calculated as described in
the text. There are multiple regions with “internal” chip-scale alignments at
the $>95\%$ confidence level. The middle row shows the average complex
ellipticity magnitude $\langle e\rangle$ (colors) and phase (lines in the
“north-up” frame) for these statistically significant bins alone. Note how
$\langle e\rangle\gtrsim 0.1$ in all significant cases suggestive of strong
alignments.

Figure 12 shows three example chips that have the highest statistical
significance ($p\lesssim 0.01$). The distribution of $\langle e\rangle$ under
the null hypothesis of random orientations can be ruled out with $\gtrsim
99\%$ confidence given the large observed $\langle e\rangle$. Coherent
alignments in the orientations of background galaxies can be seen with hints
of circular polarization patterns. Whether these alignments are due to lensing
from foreground substructure, intrinsic alignments, systematics or pure random
chance remains to be seen and will require detailed follow-up spectroscopy and
lens modeling. Appendix A shows that PSF uncertainties are unlikely to explain
these alignments. Another way to get a handle on such systematics is to search
for alignments on even larger scales which we turn to next.

Figure 12: Three example NIRCam chips in which the null hypothesis of random
orientations can be ruled out at $\gtrsim 99\%$ significance when averaging
over the combined background sample (left-most column of Figure 11). A two
arcsec rectangle is plotted at the location of each background galaxy denoting
the orientation of its major axis. The bottom three panels of each column show
the distribution of $\langle e\rangle$ under the null hypothesis of random
orientations (top), distribution of ellipticities (middle) and and
distribution of position angles (bottom). Coherent alignments in the
orientations of background galaxies can be seen.

###  5.3 Alignments in NIRCam modules $(2.2^{\prime}\times 2.2^{\prime})$

On large scales, assuming background galaxies are randomly oriented, we should
find fewer regions with statistically significant alignments unless lensing or
systematics are playing a role. Here we increase the map “pixel” size from the
scale of an individual NIRCam chip ($64^{\prime\prime}\times
64^{\prime\prime}$) to an individual NIRCam module ($2.2^{\prime}\times
2.2^{\prime}$ comprising four chips).999Note that a single NIRCam module can
sometimes fit a single bright foreground cluster used for conventional strong
and weak lensing studies (e.g., Pontoppidan et al., 2022) whereas here we are
searching for alignments over multiple modules or pointings. Figure 13
identifies a handful of NIRCam modules with statistically significant
alignments at the $>95\%$ level using our null test described in section 4.5.
A few of these have module-scale significance at the $>99\%$ level. Given that
we have 20 such modules, we expect $20\times 0.05\approx 1$ such detection by
random chance at the $95\%$ level. We have $2-3$ modules with detections for
each subsample so we pass this simple threshold. Formally, the Benjamini-
Hochberg and Bonferroni corrections lead to $p>0.05$ for all modules except
one in the $z>2$ low-mass map which remains significant at the $\sim 98\%$
level ($p=0.01882$).

Figure 14 shows the module with the most significant alignments in the
orientations of low-mass galaxies at $z>2$. Even on these larger
$2.2^{\prime}\times 2.2^{\prime}$ scales, coherent alignments can be seen.
Some of these alignments appear like linear “chains” whereas others arise from
nearly circular arrangements of background galaxies. The latter motivates
future decomposition of the shear map into E- and B-mode polarizations to
assess whether this signal is due to lensing, intrinsic alignments or
systematics.

Figure 13: Similar to Figure 11 but now for individual NIRCam modules
($2.2^{\prime}\times 2.2^{\prime}$). There are a few modules with alignments
at the $>95\%$ significance level for each subsample, with many showing
$\langle e\rangle\sim 10\%$ even on these larger scales. Figure 14: Similar to
Figure 12 but for the NIRCam module with the most statistically significant
alignments for the low-mass $z>2$ subsample (rightmost column of Figure 13).
Again coherent alignments can be seen even on these larger $2.2^{\prime}\times
2^{\prime}2$ scales, with hints of circular patterns that motivate future
$E/B$-mode polarization decomposition and weak lensing mass reconstruction.

###  5.4 Alignments in NIRCam pointings $(2.2^{\prime}\times 5.1^{\prime})$

Here we increase the map “pixel” size to even larger scales of an entire
NIRCam pointing which encompasses two $2.2^{\prime}\times 2.2^{\prime}$
modules separated by a $\sim 0.7^{\prime}$ gap. Figure 15 shows that we
continue to find $1-2$ pointings for each subsample with significant
alignments at the $>95\%$ level. Two of these are significant at the $>99\%$
level: one in the southwest corner of the high-mass map and another in the
northwest corner of the $z>2$ low-mass map. Since we have 10 pointings, we
expect $10\times 0.05\approx 0.5$ pointings with a detection by random chance
alone. We are above this simple threshold for all subsamples, but the more
formal Benjamini-Hochberg and Bonferroni methods lead to corrected $p$-values
in excess of 0.05 for all pointings. However, in addition to these
conservative statistical tests, it is still instructive to look at images of
pointings with detections.

Figure 16 shows an image of the southwest pointing with significant alignments
of high-mass galaxies at $z>1$. Large-scale alignments are clearly obvious by
eye in both modules of this pointing. Many of these galaxies also show
circular alignment patterns which reinforces the need to decompose the shear
map into $E$ and $B$ mode polarizations which can help constrain whether this
is due to lensing, intrinsic alignments or systematics. More generally, our
results motivate computing galaxy-shear and shear-shear correlation functions
in “blank” JWST deep fields like CEERS that are far away from obvious, bright
foreground clusters. Such an analysis would naturally capture the dependence
of the alignment signal as a continuous function of spatial scale and offer
constraints on cosmology and the foreground lens mass. However, this is non-
trivial and requires careful consideration of survey geometry so we leave it
for the future.

Figure 15: Same as Figure 13 but now averaging over entire NIRCam pointings
(two $2.2^{\prime}\times 2.2^{\prime}$ modules with a $\sim 0.7^{\prime}$ gap
in between). For each subsample, there are $1-2$ pointings with alignments at
$>95\%$ significance. Figure 16: Same as Figure 14 but now for the southwest
$2.2^{\prime}\times 5.1^{\prime}$ pointing in which our high-mass sources have
significant alignments at the $>99\%$ level. Large-scale alignments are
clearly present and there are hints of circular polarization patterns.

###  5.5 Alignments averaged over the entire survey

Figure 17 shows that, on the scale of the entire field ($\sim
30^{\prime}\times 6^{\prime}$), we do not find any statistically significant
alignments ($p>0.05$). This makes sense since over such a large area,
background galaxies should be randomly oriented and lensing must be weak. It
also confirms the lack of systematics in the data and our shape measurements,
at least on such large scales. Note that we also found a similar lack of net
alignments in Figure 8 for galaxy-galaxy lensing pairs with negligible shear,
which comprise most of our background galaxies distributed throughout the
entire field.

Figure 17: Same as Figure 13 but now averaging over the entire field ($\sim
30^{\prime}\times 6^{\prime}$). As expected, we do not find any statistically
significant alignments implying that the background galaxies are randomly
oriented over such large scales.

###  5.6 Shear correlations as a function of scale

Lastly, to summarize, Figure 18 shows how the average complex ellipticity
$\langle e\rangle$ and shear variance $\langle\overline{\gamma}^{2}\rangle$
depend on scale for both our background galaxies and PSF stars (see Appendix A
for more details about the latter). Recall that $\langle e\rangle$ gives the
average orientation and ellipticity of all galaxies in some region whereas
$\langle\overline{\gamma}^{2}\rangle$ quantifies the degree of correlation
between pairwise shapes (as detailed in 3.5). In this context, since
foreground stars should not be sheared (unlike high-redshift galaxies), they
serve as a null test for PSF systematics by providing a “floor” on $\langle
e\rangle$ and $\langle\overline{\gamma}^{2}\rangle$. This is complementary to
our null test for statistical significance in previous subsections.

Using our combined background galaxy sample without splitting by mass,
redshift or alignment significance, we find that $\langle e\rangle$ drops from
an average of $\sim 8\%$ at the NIRCam chip scale to $\sim 5\%$ for NIRCam
modules and $\sim 3\%$ for NIRCam pointings. If we restrict ourselves only to
chips, modules and pointings where the null hypothesis of randomly oriented
galaxies could be ruled out at $>95\%$ significance, then $\langle e\rangle$
is somewhat larger due to the alignments, as expected. The exact same
calculation for the stars gives $\langle e\rangle$ that is almost constant
with scale but drops monotonically from $\sim 2-2.5\%$ for short-wavelength
NIRCam filters (F115W, F150W, F200W) to $\sim 0.5-1.5\%$ for long-wavelength
NIRCam filters (F277W, F356W, F444W). This is several times lower than our
$\langle e\rangle$ measured for galaxies at the chip and module scales, and to
a lesser extent for pointings. This implies that the alignments we have
detected on these smaller scales cannot be attributed to percent-level PSF
systematics (see Appendix A for more justification). In contrast, on the scale
of the entire survey, the residual $\langle e\rangle\sim 0.9\%$ of galaxies is
comparable to or even exceeded by the $\langle e\rangle$ of PSF stars, which
is consistent with the alignments vanishing on large scales as they should.

The middle panel of Figure 18 similarly shows that the shear variance drops
from $\langle\overline{\gamma}^{2}\rangle\sim 4\times 10^{-3}$ at chip scales
to $10^{-3}$ for modules and $6\times 10^{-4}$ for pointings. The exact same
calculation for the PSF stars leads to $\langle\overline{\gamma}^{2}\rangle$
that is almost constant with scale and drops nearly monotonically from
$4\times 10^{-4}$ for the bluest filter (F115W) to $sim1-3\times 10^{-5}$ for
the reddest filter (F444W). Thus, pairs of galaxies have significant shear
correlations above that expected from PSF systematics on the scale of chips,
modules and pointings. In contrast, on the scale of the entire survey,
$\langle\overline{\gamma}^{2}\rangle\sim 3\times 10^{-5}$ for galaxies, which
is comparable to or even smaller than that of the stars. Thus the shear
variance on the survey scale is consistent with being zero since it is so low
as to be in the PSF systematics-dominated regime.

If our alignments were due to cosmic shear, we would expect a typical
$\langle\overline{\gamma}^{2}\rangle\sim 3\times 10^{-4}$ on chip scales that
decreases to $\sim 3\times 10^{-5}$ on the scale of the survey (based on
Figure 7 of the review by Refregier, 2003). Our measured shear variance is an
order of magnitude larger at the chip scale and several times larger at the
module and pointing scales, though consistent on the survey scale. The right
panel of Figure 18 shows that the contribution of the star-galaxy cross-
correlation is of order $\sim 10\%$ the shear variance on the chip scale,
which means that star-galaxy alignments must be weaker than galaxy-galaxy
alignments on these small scales. However, the fractional contribution of the
star-galaxy cross-correlation becomes more significant on larger scales where
the shear signal itself is lower and thus in the systematics-dominated regime.
Thus our measured shear variance is an upper limit and it implies a root-mean-
square shear $\sqrt{\langle\overline{\gamma}^{2}\rangle}\sim 5\%$ which is
still too small to explain the preferentially high ellipticity $e\sim 0.4-0.7$
of early low-mass galaxies. The actual shear is likely even smaller and future
studies that correct galaxy shapes for the local rather than global PSF may be
able to overcome systematics and measure a much weaker cosmic shear signal
expected in the data that we cannot.

Figure 18: Average complex ellipticity (left), shear variance (middle) and
star-galaxy cross-correlation (right) on the scale of individual NIRCam chips,
modules and pointings as well as the entire survey. All points denote averages
and standard errors. Solid black lines show results for galaxies in all
regions whereas dashed black lines are only based on regions where the null
hypothesis of randomly oriented galaxies could be ruled out at $>95\%$
significance. Colored lines show the exact same calculations for PSF stars in
each filter. Both $\langle e\rangle$ and $\langle\overline{\gamma}^{2}\rangle$
are above the floor set by stars on smaller scales, but drop into the
systematics-dominated regime on larger scales. The shear variance in
particular is several times larger than the typical value from cosmic shear
(gray shaded region). The fractional contribution of the star-galaxy cross-
correlation to the shear variance is sub-dominant on the chip scale but
becomes significant on larger scales where the expected shear signal is lower.
Note that we do not cross-correlate stars and galaxies in F277W because that
filter is not used. This implies that our measured shear variance is an upper
limit but even the large value on the chip scale only implies a root-mean-
square shear of $\sqrt{\langle\overline{\gamma}^{2}\rangle}\sim 5\%$ which is
much smaller than the preferentially high ellipticity $e\sim 0.4-0.7$ of early
low-mass galaxies. Future work that corrects galaxy shapes for the local
rather than global PSF may be able to measure the weak expected cosmic shear
signal in “blank” JWST deep fields.

Table 2 gives alignment summary statistics for each of the four subsamples in
all 111 regions (entire survey, 10 pointings, 20 modules, 80 chips).

Region | RA | Dec | $N^{\rm all}$ | $p^{\rm all}$ | $\langle e\rangle^{\rm all}$ | $\sigma_{\langle e\rangle}^{\rm all}$ | $\langle\phi\rangle^{\rm all}$ | $\sigma_{\langle\phi\rangle}^{\rm all}$ | $|\overline{\gamma}^{2}|$ | $\sigma_{|\overline{\gamma}^{2}|}$ |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Survey | nan | nan | 3848 | 0.351439 | 0.008596 | 0.000890 | -76.610408 | 2.856324 | 0.000028 | 0.000008 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1 | 214.984415 | 52.978435 | 452 | 0.066638 | 0.038976 | 0.003063 | -77.855550 | 2.039165 | 0.000576 | 0.000111 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_A | 214.958701 | 52.960792 | 199 | 0.720549 | 0.020503 | 0.004996 | -79.233521 | 6.838656 | 0.000463 | 0.000106 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_B | 215.010129 | 52.996078 | 253 | 0.053044 | 0.054112 | 0.003427 | -77.571261 | 1.917349 | 0.001001 | 0.000183 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_A1 | 214.940923 | 52.960792 | 36 | 0.065806 | 0.124969 | 0.009176 | 87.494525 | 2.213625 | 0.005178 | 0.001155 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_A2 | 214.958701 | 52.943014 | 49 | 0.813965 | 0.030689 | 0.009761 | -88.482837 | 83.918228 | 0.003371 | 0.000461 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_A3 | 214.958701 | 52.978570 | 53 | 0.463524 | 0.062890 | 0.013377 | -50.235209 | 5.957865 | 0.001892 | 0.000484 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
P1_A4 | 214.976479 | 52.960792 | 59 | 0.745130 | 0.037742 | 0.007163 | 40.715718 | 6.197549 | 0.001791 | 0.000253 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Table 2: Alignment summary statistics for all 111 regions (entire survey, 10
pointings, 20 modules and 80 chips). The central RA and Dec of each region are
in decimal degrees. $N$ is the number of galaxies. We give the $p$-value from
our null test, magnitude of the average complex ellipticity $\langle e\rangle$
and the associated phase angle $-90^{\circ}<\langle\phi\rangle<90^{\circ}$
defined east of north. All uncertainties are standard errors from our Monte
Carlo method. We give these values for all four subsamples denoted by the
superscript in the column name: all, high-mass, low-mass or “early” for $z>2$
low-mass galaxies. The shear variance and star-galaxy cross-correlation are
given only for the combined sample. The full table is available for download
from the journal.

##  6 Discussion

###  6.1 Evidence for cosmic shear in JWST-CEERS?

We presented evidence for large-scale alignments of high-redshift galaxies in
a “blank” JWST deep field far away from any obvious foreground cluster.
Galaxies at $z>1$ with $\log_{10}M_{*}/M_{\odot}>9$ appear to be aligned on
the scale of multiple NIRCam chips ($64^{\prime\prime}\times
64^{\prime\prime}$), modules ($2.2^{\prime}\times 2.2^{\prime}$) and pointings
($2.2^{\prime}\times 5.1^{\prime}$) throughout the JWST-CEERS survey
footprint. When averaging over background galaxies in these regions, the norm
of the complex ellipticity vector is often $\langle e\rangle\gtrsim 10\%$,
which is at least an order of magnitude above the conventional weak lensing
regime. Many of these chips, modules and pointings reveal circular alignment
patterns that are suggestive of non-random effects and motivate follow-up
$E/B$-mode polarization decomposition. In addition, our tangentially-aligned
galaxy-galaxy lensing candidates, including one possible strong lens, tend to
cluster in the southwest (Figure 9). This clustering of galaxy-galaxy lensing
candidates on scales of $\sim 0.1$ deg as well as the detection of alignments
in multiple regions on scales up to $2.2^{\prime}\times 5.1^{\prime}$ equates
to several Mpc at $z\sim 0.75$ where we found an overdensity of massive
foreground lens galaxies (Figure 1). It is plausible that a cluster, proto-
cluster or filament at $z\sim 0.75$ may be responsible for the lensing in
which case detailed follow-up modeling may help us learn something about the
intrinsic 3D geometry and lensing mass profile of this overdensity (e.g.,
Kovner, 1987; Keeton et al., 1997).

On the other hand, we cannot rule out intrinsic alignments in the background
galaxies (Kirk et al., 2015; Mandelbaum, 2018; Lamman et al., 2024). Pandya et
al. (2019) showed that if early low-mass galaxies are intrinsically prolate,
then they are expected to show strong intrinsic alignments with each other as
they form along cosmic web filaments. Since our analysis is not done in bins
of source redshift, we are presumably averaging over this and therefore our
signal should not be driven by intrinsic alignments. Pandya et al. (2019) did
not detect the expected intrinsic alignments and claimed that it was due to a
lack of spectroscopically-confirmed background pairs in HST-CANDELS. However,
they also did not use the complex ellipicity formalism employed here so the
issue is worth revisiting in light of our results. The fact that we are
detecting coherent alignments on the scale of entire NIRCam pointings
($2.2^{\prime}\times 5.1^{\prime}$) and finding $\langle e\rangle\gtrsim 10\%$
in multiple areas spread across the entire $\sim 0.3$ deg extent of the survey
(several Mpc at $z\sim 0.75$) does point towards lensing. But because we lack
spectroscopic redshifts for most of our background galaxies, we cannot
precisely identify high-redshift pairs to robustly constrain intrinsic
alignments.

Of course, there could be other systematics at play. It is unlikely that we
are contaminated by satellite-central pair alignments (Hirata et al., 2004)
because moderate/strong galaxy-galaxy shear requires maximizing the angular
diameter distance between the source and lens which implies very different
photometric redshift distributions. Appendix A shows that our PSF systematics
are at the few percent-level which, while inadequate for traditional weak
lensing, are an order of magnitude lower than our $\langle e\rangle\gtrsim
10\%$ from averaging over the orientations of background galaxies on the chip-
scale. Thus PSF uncertainties are unlikely to explain our alignments in terms
of $\langle e\rangle$ but are prohibitively large for measuring the precise
amplitude of the cosmic shear variance signal
$\langle\overline{\gamma}^{2}\rangle$ (Figure 18). Astrometric misalignments
could cause distortions that produce spurious shear signals (Finner et al.,
2023b) but the CEERS imaging we use is already tied to Gaia (Bagley et al.,
2023) and we have visually inspected all of our galaxies so this is unlikely.
Future decomposition of the shear field into $E$ and $B$ mode polarizations
can help constrain some of these systematics including signatures of intrinsic
alignments.

We note that JWST-CEERS is only the latest in a long line of multi-wavelength
observations of EGS (Groth et al., 1994) whose early HST-WFPC2 observations
enabled the first detection of cosmic shear from space in this very same field
(Rhodes et al., 2000, 2001). Those authors found mean shear values similar to
ours of $\sim 10\%$ using $\sim 50$ background galaxies in many $1.27\times
1.27$ arcmin2 regions along the strip (see Figure 1 of Rhodes et al., 2001).
After carefully subtracting off systematics, they estimated the shear variance
to be $\langle\overline{\gamma}^{2}\rangle\sim 0.018^{2}\sim 3\times 10^{-4}$
on the scale of a single HST-WFPC2 chip ($1.27^{\prime}\times 1.27^{\prime}$).
This is an order of magnitude lower than our
$\langle\overline{\gamma}^{2}\rangle\sim 4\times 10^{-3}$ on the slightly
smaller scale of a NIRCam chip ($64^{\prime\prime}\times 64^{\prime\prime}$).
The discrepancy is likely due to PSF systematics on our end (Figure 18) which
is why we consider our measurement of shear variance to be an upper limit
rather than a detection. However, differences in source selection and rest-
frame optical versus rest-frame ultraviolet morphology (HST-WFPC2 probed the
latter) may also play a role. Future efforts to correct for the complicated
spatially-varying PSF of JWST should enable precise measurements of the
amplitude of the cosmic shear signal in multiple “blank” JWST deep fields.
These can then be averaged to mitigate cosmic variance and help constrain the
cosmological parameters $\sigma_{8}$ and $\Omega_{\rm M}$ (as has been done
from the ground and with HST; Kaiser et al., 2000; Bacon et al., 2000; Wittman
et al., 2000; Rhodes et al., 2001; Pirzkal et al., 2001; Miralles et al.,
2002; Hämmerle et al., 2002; Rhodes et al., 2004; Schrabback et al., 2007).
Indeed, the shallower but $\sim 10\times$ larger COSMOS-Web survey is designed
in part to enable this kind of measurement with JWST (Casey et al., 2023),
building off of the long history of cosmic shear studies in the HST-COSMOS
field (Rhodes et al., 2007; Leauthaud et al., 2007; Massey et al., 2007b).

EGS has benefited greatly from dedicated panchromatic imaging and
spectroscopic redshift surveys such as AEGIS, DEEP2 and DEEP3 (Davis et al.,
2003; Faber et al., 2003; Davis et al., 2007; Cooper et al., 2007, 2012;
Newman et al., 2013). This wealth of data led previous groups to identify a
number of groups and clusters in EGS between $z\sim 0.4-1.4$, some of which
are located at $z\sim 0.75$ consistent with our claimed overdensity of lenses
there (e.g., see Figure 6 of Gerke et al., 2012). The galaxy color-magnitude
bimodality and color-density relations are also observed for EGS galaxies at
$z\sim 0.75$ which suggests the presence of rich groups (Faber et al., 2007;
Cooper et al., 2007, and see also later work by Pandya et al. 2017a using HST-
CANDELS). The abundance and clustering of X-ray-selected active galactic
nuclei also point towards a number of massive groups in EGS at $z\sim 0.7-1.4$
(Coil et al., 2009). We note that our most massive lens has a spectroscopic
redshift of $z\sim 0.78$ and although its $\sigma_{*}\sim 580$ km/s is very
high compared to that of the most massive local galaxies (Ma et al., 2014;
Pandya et al., 2017b; Greene et al., 2019), it seems consistent with the tail
of high $\sigma_{*}$ LRGs out to $z\sim 0.7$ (Thomas et al., 2013). This
massive lens is also associated with our only strong lensing candidate (Figure
5). Spectroscopic constraints on the stellar kinematics and environment of
this massive galaxy combined with secure redshifts for background sources
would enable more sophisticated source-lens reconstruction. Finally, Moustakas
et al. (2007) conducted a visual search for strong lenses in EGS using $\sim
650$ arcmin2 of HST-ACS imaging. They found three (an Einstein Cross, one with
two pairs of arcs, and one with a single pair of arcs) but unfortunately those
are all outside of the JWST-CEERS footprint which covers only $\sim 100$
arcmin2 of EGS. However, they also found that each of these lenses required
$\sim 10\%$ external shear which supports our claim that there may be
significant cosmic shear throughout CEERS leading to $\langle e\rangle\sim
0.1$.

Our study suggests that “blank” JWST deep fields may provide complementary
constraints on lensing and intrinsic alignments beyond what is traditionally
available from massive cluster surveys, provided that the systematics are
better understood. For example, the high resolution of JWST imaging may allow
one to study higher order distortions due to lensing such as flexion (which,
to quote Schneider & Er 2008, happens when “weak lensing goes bananas”). The
prospect of searching for overdensities in “blank” fields is tantalizing
because one of the strengths of lensing is that it is sensitive to mass
concentrations even if it is “dark.” JWST will also enable the selection of
lower mass, higher redshift sources for weak lensing which, while improving
number statistics, may introduce new biases due to the unknown intrinsic
shapes and intrinsic alignments of the early galaxy population. On the galaxy-
galaxy lensing side, Holloway et al. (2023) used mock data to show that JWST
deep fields should expect to have $\sim 25-65$ strongly lensed systems. One
new Einstein ring has already been detected in COSMOS-Web (van Dokkum et al.,
2024; Mercier et al., 2023). We verified that if we had a source-lens pair in
CEERS with parameters similar to the COSMOS-Web Einstein ring, we would have
identified it. Assuming rings/arcs/clumps were optimally deblended, this
implies there are no Einstein rings in CEERS.101010We repeated our search
allowing for massive lenses at any redshift and background sources down to
$\log_{10}M_{*}/M_{\odot}=8$ but did not find any convincing strong lensing
candidates in JWST-CEERS. Given the small sizes and faint magnitudes of many
of our tangentially-aligned source-lens pairs (Figure 7), follow-up
observations may be challenging but required to enable confirmation and more
detailed lensing modeling.

Looking to the near future, the Roman Space Telescope, with its $\sim 0.28$
deg2 field of view and $\sim 0.1$ arcsec resolution, has the potential to
dramatically transform our understanding of lensing, intrinsic alignments and
origin of low-mass elongated galaxies. A single pointing with Roman would
image a region larger than the entirety of the existing HST ACS/WFC3 footprint
in EGS, which required tens of tiles over only $\sim 0.05-0.2$ deg2 (Stefanon
et al., 2017). A “Roman Deep Field” (Koekemoer et al., 2019; Drakos et al.,
2022) that can reach $\sim 27$ AB mag at $\sim 1.5\mu$m would detect and
resolve elongated $\log M_{*}/M_{\odot}>9$ galaxies out to $z\sim 3$ (Pandya
et al., 2024). This would uniquely enable the calculation of galaxy-shear and
shear-shear correlation functions over much larger areas than possible with
HST or JWST. Multiple such Roman deep fields spread across the sky would help
mitigate the cosmic variance that has plagued the field for decades (Madau &
Dickinson, 2014) and allow us to test whether the alignments are indeed due to
lensing since the foreground mass distribution should vary across the sky. We
caution that the current plan for the Roman High Latitude Wide Area Survey
over $\sim 2000$ deg2 is only designed to reach $26.5$ AB mag for point
sources (Troxel et al., 2021; Wang et al., 2022; Montes et al., 2023). A rule
of thumb is that the completeness limit for extended high-redshift galaxies is
$\sim 2$ magnitudes fainter than the point source limit (Pandya et al., 2024)
so a deeper imaging component is necessary for our science case. Note that
Euclid, with its sensitivity and $\sim 0.2$ arcsec resolution, is not expected
to be complete to such low-mass galaxies beyond $z\sim 0.5$ even in its Deep
Surveys spanning $\sim 50$ deg2 (Euclid Collaboration et al., 2022).

###  6.2 Why are early galaxies preferentially elongated?

Despite the evidence for galaxy-galaxy lensing and large-scale alignments with
$\langle e\rangle\sim 10\%$, the fact that most of our background sources have
much higher ellipticities of $e\sim 0.4-0.7$ (top panel of Figure 3) implies
that gravitational lensing cannot be primarily responsible for their
elongation. Thus, we are led back to the original question that motivated this
study: why are early low-mass galaxies preferentially elongated as observed by
HST (Cowie et al., 1995; Elmegreen et al., 2005; Ravindranath et al., 2006;
van der Wel et al., 2014; Zhang et al., 2019) and JWST (Pandya et al., 2024)?

The remaining explanations for this puzzle range from the mundane to the
exotic. There may still be a surface brightness detection bias against
rounder, face-on disks (Loeb, 2024). Pandya et al. (2024) claimed that JWST-
CEERS is sufficiently complete but deeper surveys have not yet been checked.
Relatedly, we may be seeing elongated star-forming proto-bars at the centers
of extremely low surface brightness disks. Alternatively, the elongated shapes
reflect tidal debris from ongoing mergers. Finally, early low-mass galaxies
may not start out as disks as commonly assumed but rather as triaxial/prolate
structures that reflect the cosmic web filaments in which they form (Ceverino
et al., 2015; Tomassetti et al., 2016; Pandya et al., 2019). Pandya et al.
(2024) emphasized that the stellar masses of early elongated galaxies at
$z\sim 2-3$ are consistent with them being Milky Way progenitors, broadly
defined, implying that our own Galaxy may have been elongated in its past.
There are hints of this early elongation from Galactic archaeology (Naidu et
al., 2021), and the recent discovery of a strongly lensed, elongated, clumpy
Milky Way-like progenitor at $z\sim 8.3$ evokes connections to chain galaxies
(Mowla et al., 2024, see also Zonoozi et al. 2019).

It is imperative that we try to rule out the prolate/triaxial interpretation
because so much of our current understanding of galaxy formation depends on
them starting out as disks due to tidal torques and angular momentum
conservation (Mo et al., 1998). Although prolate galaxies can form in
$\Lambda$CDM, this has only been seen in some cosmological zoom-in simulations
(Ceverino et al., 2015; Tomassetti et al., 2016) and there is no consensus on
subgrid models for uncertain baryonic physics (e.g., Somerville & Davé, 2015).
What we really want to know is the large-scale clustering, orientation,
dynamical stability and cosmological evolution of prolate galaxies but large-
volume simulations typically do not have the resolution required to robustly
predict the internal dynamics and therefore 3D shapes of early low-mass
galaxies (e.g., Pillepich et al., 2019). All of this invites an exploration of
more exotic scenarios such as alternatives to cold dark matter. For example,
it has already been shown that warm and axionic (fuzzy) dark matter naturally
leads to more elongated halos and therefore possibly more elongated galaxies
at early times as observed (Mocz et al., 2020; Dome et al., 2023).

Pandya et al. (2024) claimed that stellar kinematics will be required to
definitively distinguish between these interpretations. However, our analysis
suggests a new test: if galaxy-shear and shear-shear correlation functions can
be measured over large deep fields for $\log_{10}(M_{*}/M_{\odot})\sim 9$
galaxies out to $z\sim 2$, then we will be able to constrain whether these
early elongated galaxies show strong intrinsic alignments and trace cosmic web
filaments at high-redshift as predicted by Pandya et al. (2019). The Euclid
Deep Surveys are not expected to detect and resolve such low-mass galaxies
beyond $z\sim 0.5$ (Euclid Collaboration et al., 2022), but the Roman Space
Telescope may be able to do so. We therefore strongly recommend that a Roman
Deep Field (Koekemoer et al., 2019) be planned to enable this potentially
transformative science: are early elongated galaxies reliable probes of
lensing, do they represent an unrecognized source of significant additive
bias, and can they tell us anything fundamentally new about cosmology?

##  7 Summary

We used data from JWST-CEERS (Finkelstein et al., 2023, McGrath et al., in
prep.) to explore an idea not previously considered in the literature: are
early low-mass galaxies preferentially elongated because of gravitational
lensing? First, we investigated galaxy-galaxy lensing for which we selected 76
massive foreground lens galaxies at $z<1$ and 3848 background sources at $z>1$
with $\log_{10}M_{*}/M_{\odot}>9$. Above this mass limit, CEERS has, on
average, $\sim 50-100$ background galaxies per arcmin2 over its $\sim 100$
arcmin2 footprint, making it a good testbed for lensing studies in “blank”
JWST deep fields. For each source, we assigned the nearest on-sky massive lens
and estimated the shear $\gamma$ assuming a singular isothermal sphere lens
model as well as the tangential alignment angle. Second, we looked for
alignments on larger scales by averaging over the orientations of background
galaxies on the scale of individual NIRcam chips ($64^{\prime\prime}\times
64^{\prime\prime}$), modules ($2.2^{\prime}\times 2.2^{\prime}$) and pointings
($2.2^{\prime}\times 5.1^{\prime}$) across the survey. Our results are
summarized as follows.

  1. 1.

Galaxy-galaxy lensing cannot be the primary driver for the elongated
appearance of most early low-mass galaxies because our typical predicted shear
is of order only $\gamma\approx 1\%$. This is much lower than the typical
ellipticity $e\sim 0.4-0.7$ of high-redshift galaxies (Figure 3). However,
there is a broad tail towards larger shears in excess of $10\%$. We use non-
parametric quantile regression with Bayesian Additive Regression Trees to show
that there is a hint of an excess of tangentially-aligned source-lens pairs
with such high shear (Figure 4). We identify 17 tangentially-aligned lensing
candidates, only one of which may truly be in the strong lensing regime
(Figure 5 and Table 1). These galaxy-galaxy lensing candidates are clustered
in the southwest region of the survey (Figure 9) and show evidence that their
orientations are correlated (Figure 8).

  2. 2.

We find evidence for coherent large-scale alignments when averaging over the
orientations of background galaxies in multiple NIRCam chips
($64^{\prime\prime}\times 64^{\prime\prime}$), modules ($2.2^{\prime}\times
2.2^{\prime}$) and pointings ($2.2^{\prime}\times 5.1^{\prime}$). We can rule
out the null hypothesis of random orientations at the $\gtrsim 99\%$
confidence level for multiple individual regions (Figures 11, 13 and 15, and
Table 2). The number of such regions is small on the scale of the survey so
formally the detections may be due to random chance, but the imaging clearly
reveals non-random alignment patterns (Figures 12, 14 and 16). On the chip
scale, the average complex ellipticity $\langle e\rangle\gtrsim 10\%$ and
shear variance $\langle\overline{\gamma}^{2}\rangle\sim 4\times 10^{-3}$ are
an order of magnitude above the conventional weak lensing regime (Figure 18).
We consider these to be upper limits due to PSF uncertainties (Appendix A),
intrinsic alignments, cosmic variance and other lensing systematics. The
maximum implied “cosmic shear” is still only
$\sqrt{\langle\overline{\gamma}^{2}\rangle}\sim 5\%$, much lower than the
ellipticity of early low-mass galaxies.

  3. 3.

We speculate that the large-scale alignments may be due to lensing from a
foreground protocluster or filament spanning several Mpc at $z\sim 0.75$ where
we found an overabundance of massive lens galaxies (Figure 1). More detailed
follow-up weak lensing analyses in CEERS and other “blank” JWST deep fields
are needed to confirm whether the alignments are indeed due to lensing by
large-scale structure, intrinsic alignments or other lensing systematics. We
also strongly recommend the planning of a deep field with the forthcoming
Roman Space Telescope whose sensitivity, resolution and $0.28$ deg2 field of
view will uniquely enable the calculation of galaxy-shear and shear-shear
correlation functions for elongated galaxies with
$\log_{10}(M_{*}/M_{\odot})\sim 9$ out to $z\sim 2$. This will provide a novel
way to constrain the intrinsic shapes of early elongated galaxies which, if
they are forming along filaments, are expected to show strong intrinsic
alignments (Pandya et al., 2019). This would also help calibrate any bias from
this population for weak lensing analyses.

This paper is the result of an initial scientific disagreement between VP and
AL about whether “galaxies are going bananas” (Pandya et al., 2024; Loeb,
2024). Thanks to an invitation from AL, VP gave a colloquium at Harvard’s ITC
and that visit inspired conversations about the enigmatic nature of early
elongated galaxies and novel ways to follow them up. We thank Haowen Zhang,
Raphael Gavazzi and Greg Bryan for helpful suggestions. VP thanks Jo Bovy for
writing the interactive online “Dynamics and Astrophysics of Galaxies”
textbook at https://galaxiesbook.org/. VP also thanks Joel Primack, Avishai
Dekel, Sandy Faber and David Koo for getting him interested in early elongated
galaxies. Support for VP was provided by NASA through the NASA Hubble
Fellowship grant HST-HF2-51489 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555.

## References

  * Bacon et al. (2000) Bacon, D. J., Refregier, A. R., & Ellis, R. S. 2000, MNRAS, 318, 625, doi: 10.1046/j.1365-8711.2000.03851.x
  * Bagley et al. (2023) Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJ, 946, L12, doi: 10.3847/2041-8213/acbb08
  * Barkana & Loeb (2000) Barkana, R., & Loeb, A. 2000, ApJ, 531, 613, doi: 10.1086/308503
  * Barro et al. (2023) Barro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2305.14418, doi: 10.48550/arXiv.2305.14418
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bolton et al. (2008) Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, ApJ, 682, 964, doi: 10.1086/589327
  * Bradley et al. (2024) Bradley, L. D., Adamo, A., Vanzella, E., et al. 2024, arXiv e-prints, arXiv:2404.10770, doi: 10.48550/arXiv.2404.10770
  * Brainerd et al. (1996) Brainerd, T. G., Blandford, R. D., & Smail, I. 1996, ApJ, 466, 623, doi: 10.1086/177537
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Bridle et al. (2010) Bridle, S., Balan, S. T., Bethge, M., et al. 2010, MNRAS, 405, 2044, doi: 10.1111/j.1365-2966.2010.16598.x
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597, doi: 10.1146/annurev-astro-082214-122432
  * Casey et al. (2023) Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31, doi: 10.3847/1538-4357/acc2bc
  * Ceverino et al. (2015) Ceverino, D., Primack, J., & Dekel, A. 2015, MNRAS, 453, 408, doi: 10.1093/mnras/stv1603
  * Cha et al. (2024) Cha, S., HyeongHan, K., Scofield, Z. P., Joo, H., & Jee, M. J. 2024, ApJ, 961, 186, doi: 10.3847/1538-4357/ad0cbf
  * Chang et al. (2013) Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013, ApJ, 773, 149, doi: 10.1088/0004-637X/773/2/149
  * Chipman et al. (2010) Chipman, H. A., George, E. I., & McCulloch, R. E. 2010, The Annals of Applied Statistics, 4, doi: 10.1214/09-aoas285
  * Coil et al. (2009) Coil, A. L., Georgakakis, A., Newman, J. A., et al. 2009, ApJ, 701, 1484, doi: 10.1088/0004-637X/701/2/1484
  * Conselice (2014) Conselice, C. J. 2014, ARA&A, 52, 291, doi: 10.1146/annurev-astro-081913-040037
  * Cooper et al. (2007) Cooper, M. C., Newman, J. A., Coil, A. L., et al. 2007, MNRAS, 376, 1445, doi: 10.1111/j.1365-2966.2007.11534.x
  * Cooper et al. (2012) Cooper, M. C., Griffith, R. L., Newman, J. A., et al. 2012, MNRAS, 419, 3018, doi: 10.1111/j.1365-2966.2011.19938.x
  * Cowie et al. (1995) Cowie, L. L., Hu, E. M., & Songaila, A. 1995, AJ, 110, 1576, doi: 10.1086/117631
  * Dalcanton & Shectman (1996) Dalcanton, J. J., & Shectman, S. A. 1996, ApJ, 465, L9, doi: 10.1086/310124
  * Davis et al. (2003) Davis, M., Faber, S. M., Newman, J., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4834, Discoveries and Research Prospects from 6- to 10-Meter-Class Telescopes II, ed. P. Guhathakurta, 161–172, doi: 10.1117/12.457897
  * Davis et al. (2007) Davis, M., Guhathakurta, P., Konidaris, N. P., et al. 2007, ApJ, 660, L1, doi: 10.1086/517931
  * Dome et al. (2023) Dome, T., Fialkov, A., Mocz, P., et al. 2023, MNRAS, 519, 4183, doi: 10.1093/mnras/stac3766
  * Drakos et al. (2022) Drakos, N. E., Villasenor, B., Robertson, B. E., et al. 2022, ApJ, 926, 194, doi: 10.3847/1538-4357/ac46fb
  * Elmegreen et al. (2005) Elmegreen, D. M., Elmegreen, B. G., Rubin, D. S., & Schaffer, M. A. 2005, ApJ, 631, 85, doi: 10.1086/432502
  * Euclid Collaboration et al. (2022) Euclid Collaboration, Bretonnière, H., Huertas-Company, M., et al. 2022, A&A, 657, A90, doi: 10.1051/0004-6361/202141393
  * Faber et al. (2003) Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1657–1669, doi: 10.1117/12.460346
  * Faber et al. (2007) Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265, doi: 10.1086/519294
  * Feigelson & Babu (2012) Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods for Astronomy, doi: 10.48550/arXiv.1205.2064
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Finner et al. (2023a) Finner, K., Faisst, A., Chary, R.-R., & Jee, M. J. 2023a, ApJ, 953, 102, doi: 10.3847/1538-4357/ace1e6
  * Finner et al. (2023b) Finner, K., Lee, B., Chary, R.-R., et al. 2023b, ApJ, 958, 33, doi: 10.3847/1538-4357/acfafd
  * Frye et al. (2023) Frye, B. L., Pascale, M., Foo, N., et al. 2023, ApJ, 952, 81, doi: 10.3847/1538-4357/acd929
  * Gerke et al. (2012) Gerke, B. F., Newman, J. A., Davis, M., et al. 2012, ApJ, 751, 50, doi: 10.1088/0004-637X/751/1/50
  * Greene et al. (2019) Greene, J. E., Veale, M., Ma, C.-P., et al. 2019, ApJ, 874, 66, doi: 10.3847/1538-4357/ab01e3
  * Groth et al. (1994) Groth, E. J., Kristian, J. A., Lynds, R., et al. 1994, in American Astronomical Society Meeting Abstracts, Vol. 185, American Astronomical Society Meeting Abstracts, 53.09
  * Hämmerle et al. (2002) Hämmerle, H., Miralles, J. M., Schneider, P., et al. 2002, A&A, 385, 743, doi: 10.1051/0004-6361:20020195
  * Harvey & Massey (2024) Harvey, D. R., & Massey, R. 2024, MNRAS, 529, 802, doi: 10.1093/mnras/stae370
  * Heymans et al. (2005) Heymans, C., Brown, M. L., Barden, M., et al. 2005, MNRAS, 361, 160, doi: 10.1111/j.1365-2966.2005.09152.x
  * Heymans et al. (2006) Heymans, C., Van Waerbeke, L., Bacon, D., et al. 2006, MNRAS, 368, 1323, doi: 10.1111/j.1365-2966.2006.10198.x
  * Hill et al. (2020) Hill, J., Linero, A., & Murray, J. 2020, Annual Review of Statistics and Its Application, 7, 251, doi: https://doi.org/10.1146/annurev-statistics-031219-041110
  * Hirata et al. (2004) Hirata, C. M., Mandelbaum, R., Seljak, U., et al. 2004, MNRAS, 353, 529, doi: 10.1111/j.1365-2966.2004.08090.x
  * Holloway et al. (2023) Holloway, P., Verma, A., Marshall, P. J., More, A., & Tecza, M. 2023, MNRAS, 525, 2341, doi: 10.1093/mnras/stad2371
  * Jee et al. (2007) Jee, M. J., Blakeslee, J. P., Sirianni, M., et al. 2007, PASP, 119, 1403, doi: 10.1086/524849
  * Kaiser et al. (1995) Kaiser, N., Squires, G., & Broadhurst, T. 1995, ApJ, 449, 460, doi: 10.1086/176071
  * Kaiser et al. (2000) Kaiser, N., Wilson, G., & Luppino, G. A. 2000, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0003338
  * Kartaltepe et al. (2023) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15, doi: 10.3847/2041-8213/acad01
  * Keeton et al. (1997) Keeton, C. R., Kochanek, C. S., & Seljak, U. 1997, ApJ, 482, 604, doi: 10.1086/304172
  * Kilbinger (2015) Kilbinger, M. 2015, Reports on Progress in Physics, 78, 086901, doi: 10.1088/0034-4885/78/8/086901
  * Kirk et al. (2015) Kirk, D., Brown, M. L., Hoekstra, H., et al. 2015, Space Sci. Rev., 193, 139, doi: 10.1007/s11214-015-0213-4
  * Koekemoer et al. (2019) Koekemoer, A., Foley, R. J., Spergel, D. N., et al. 2019, BAAS, 51, 550, doi: 10.48550/arXiv.1903.06154
  * Kovner (1987) Kovner, I. 1987, ApJ, 312, 22, doi: 10.1086/164846
  * Kravtsov (2013) Kravtsov, A. V. 2013, ApJ, 764, L31, doi: 10.1088/2041-8205/764/2/L31
  * Lamman et al. (2024) Lamman, C., Tsaprazi, E., Shi, J., et al. 2024, The Open Journal of Astrophysics, 7, 14, doi: 10.21105/astro.2309.08605
  * Leauthaud et al. (2007) Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219, doi: 10.1086/516598
  * Loeb (2024) Loeb, A. 2024, arXiv e-prints, arXiv:2401.03592, doi: 10.48550/arXiv.2401.03592
  * Ma et al. (2014) Ma, C.-P., Greene, J. E., McConnell, N., et al. 2014, ApJ, 795, 158, doi: 10.1088/0004-637X/795/2/158
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Mandelbaum (2018) Mandelbaum, R. 2018, ARA&A, 56, 393, doi: 10.1146/annurev-astro-081817-051928
  * Maoli et al. (2001) Maoli, R., Van Waerbeke, L., Mellier, Y., et al. 2001, A&A, 368, 766, doi: 10.1051/0004-6361:20010058
  * Martin et al. (2021) Martin, O. A., Kumar, R., & Lao, J. 2021, Bayesian Modeling and Computation in Python (Boca Raton)
  * Mason et al. (2015) Mason, C. A., Treu, T., Schmidt, K. B., et al. 2015, ApJ, 805, 79, doi: 10.1088/0004-637X/805/1/79
  * Massey et al. (2007a) Massey, R., Heymans, C., Bergé, J., et al. 2007a, MNRAS, 376, 13, doi: 10.1111/j.1365-2966.2006.11315.x
  * Massey et al. (2007b) Massey, R., Rhodes, J., Leauthaud, A., et al. 2007b, ApJS, 172, 239, doi: 10.1086/516599
  * Mercier et al. (2023) Mercier, W., Shuntov, M., Gavazzi, R., et al. 2023, arXiv e-prints, arXiv:2309.15986, doi: 10.48550/arXiv.2309.15986
  * Miralles et al. (2002) Miralles, J. M., Erben, T., Hämmerle, H., et al. 2002, A&A, 388, 68, doi: 10.1051/0004-6361:20020426
  * Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319, doi: 10.1046/j.1365-8711.1998.01227.x
  * Mocz et al. (2020) Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2020, MNRAS, 494, 2027, doi: 10.1093/mnras/staa738
  * Montes et al. (2023) Montes, M., Annibali, F., Bellazzini, M., et al. 2023, arXiv e-prints, arXiv:2306.09414, doi: 10.48550/arXiv.2306.09414
  * Moustakas et al. (2007) Moustakas, L. A., Marshall, P., Newman, J. A., et al. 2007, ApJ, 660, L31, doi: 10.1086/517930
  * Mowla et al. (2022) Mowla, L., Iyer, K. G., Desprez, G., et al. 2022, ApJ, 937, L35, doi: 10.3847/2041-8213/ac90ca
  * Mowla et al. (2024) Mowla, L., Iyer, K., Asada, Y., et al. 2024, arXiv e-prints, arXiv:2402.08696, doi: 10.48550/arXiv.2402.08696
  * Naidu et al. (2021) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2021, ApJ, 923, 92, doi: 10.3847/1538-4357/ac2d2d
  * Narayan & Bartelmann (1996) Narayan, R., & Bartelmann, M. 1996, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/9606001
  * Newman et al. (2013) Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5, doi: 10.1088/0067-0049/208/1/5
  * Oguri et al. (2012) Oguri, M., Bayliss, M. B., Dahle, H., et al. 2012, MNRAS, 420, 3213, doi: 10.1111/j.1365-2966.2011.20248.x
  * Oriol et al. (2023) Oriol, A.-P., Virgile, A., Colin, C., et al. 2023, PeerJ Computer Science, 9, e1516, doi: 10.7717/peerj-cs.1516
  * Pandya et al. (2017a) Pandya, V., Brennan, R., Somerville, R. S., et al. 2017a, MNRAS, 472, 2054, doi: 10.1093/mnras/stx2027
  * Pandya et al. (2017b) Pandya, V., Greene, J. E., Ma, C.-P., et al. 2017b, ApJ, 837, 40, doi: 10.3847/1538-4357/aa5ebc
  * Pandya et al. (2019) Pandya, V., Primack, J., Behroozi, P., et al. 2019, MNRAS, 488, 5580, doi: 10.1093/mnras/stz2129
  * Pandya et al. (2024) Pandya, V., Zhang, H., Huertas-Company, M., et al. 2024, ApJ, 963, 54, doi: 10.3847/1538-4357/ad1a13
  * Pascale et al. (2022) Pascale, M., Frye, B. L., Diego, J., et al. 2022, ApJ, 938, L6, doi: 10.3847/2041-8213/ac9316
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196, doi: 10.1093/mnras/stz2338
  * Pirzkal et al. (2001) Pirzkal, N., Collodel, L., Erben, T., et al. 2001, A&A, 375, 351, doi: 10.1051/0004-6361:20010836
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14, doi: 10.3847/2041-8213/ac8a4e
  * Quiroga et al. (2023) Quiroga, M., Garay, P. G., Alonso, J. M., Loyola, J. M., & Martin, O. A. 2023, Bayesian additive regression trees for probabilistic programming.  https://arxiv.org/abs/2206.03619
  * Ravindranath et al. (2006) Ravindranath, S., Giavalisco, M., Ferguson, H. C., et al. 2006, ApJ, 652, 963, doi: 10.1086/507016
  * Refregier (2003) Refregier, A. 2003, ARA&A, 41, 645, doi: 10.1146/annurev.astro.41.111302.102207
  * Rhodes et al. (2004) Rhodes, J., Refregier, A., Collins, N. R., et al. 2004, ApJ, 605, 29, doi: 10.1086/382181
  * Rhodes et al. (2000) Rhodes, J., Refregier, A., & Groth, E. J. 2000, ApJ, 536, 79, doi: 10.1086/308902
  * Rhodes et al. (2001) —. 2001, ApJ, 552, L85, doi: 10.1086/320336
  * Rhodes et al. (2007) Rhodes, J. D., Massey, R. J., Albert, J., et al. 2007, ApJS, 172, 203, doi: 10.1086/516592
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Robertson et al. (2023) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, ApJ, 942, L42, doi: 10.3847/2041-8213/aca086
  * Schneider (2005) Schneider, P. 2005, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0509252
  * Schneider (2015) —. 2015, Extragalactic Astronomy and Cosmology: An Introduction, doi: 10.1007/978-3-642-54083-7
  * Schneider et al. (1992) Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses, doi: 10.1007/978-3-662-03758-4
  * Schneider & Er (2008) Schneider, P., & Er, X. 2008, A&A, 485, 363, doi: 10.1051/0004-6361:20078631
  * Schrabback et al. (2007) Schrabback, T., Erben, T., Simon, P., et al. 2007, A&A, 468, 823, doi: 10.1051/0004-6361:20065898
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Somerville et al. (2018) Somerville, R. S., Behroozi, P., Pandya, V., et al. 2018, MNRAS, 473, 2714, doi: 10.1093/mnras/stx2040
  * Stefanon et al. (2017) Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32, doi: 10.3847/1538-4365/aa66cb
  * Thomas et al. (2013) Thomas, D., Steele, O., Maraston, C., et al. 2013, MNRAS, 431, 1383, doi: 10.1093/mnras/stt261
  * Tomassetti et al. (2016) Tomassetti, M., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 458, 4477, doi: 10.1093/mnras/stw606
  * Treu (2010) Treu, T. 2010, ARA&A, 48, 87, doi: 10.1146/annurev-astro-081309-130924
  * Troxel et al. (2021) Troxel, M. A., Long, H., Hirata, C. M., et al. 2021, MNRAS, 501, 2044, doi: 10.1093/mnras/staa3658
  * Turner et al. (1984) Turner, E. L., Ostriker, J. P., & Gott, J. R., I. 1984, ApJ, 284, 1, doi: 10.1086/162379
  * Tyson et al. (1990) Tyson, J. A., Valdes, F., & Wenk, R. A. 1990, ApJ, 349, L1, doi: 10.1086/185636
  * van den Bergh et al. (1996) van den Bergh, S., Abraham, R. G., Ellis, R. S., et al. 1996, AJ, 112, 359, doi: 10.1086/118020
  * van der Wel et al. (2014) van der Wel, A., Chang, Y.-Y., Bell, E. F., et al. 2014, ApJ, 792, L6, doi: 10.1088/2041-8205/792/1/L6
  * van Dokkum et al. (2024) van Dokkum, P., Brammer, G., Wang, B., Leja, J., & Conroy, C. 2024, Nature Astronomy, 8, 119, doi: 10.1038/s41550-023-02103-9
  * Vega-Ferrero et al. (2024) Vega-Ferrero, J., Huertas-Company, M., Costantin, L., et al. 2024, ApJ, 961, 51, doi: 10.3847/1538-4357/ad05bb
  * Wang et al. (2022) Wang, Y., Zhai, Z., Alavi, A., et al. 2022, ApJ, 928, 1, doi: 10.3847/1538-4357/ac4973
  * Wittman et al. (2000) Wittman, D. M., Tyson, J. A., Kirkman, D., Dell’Antonio, I., & Bernstein, G. 2000, Nature, 405, 143, doi: 10.1038/35012001
  * Wyithe et al. (2011) Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. 2011, Nature, 469, 181, doi: 10.1038/nature09619
  * Yu & Moyeed (2001) Yu, K., & Moyeed, R. A. 2001, Statistics & Probability Letters, 54, 437, doi: https://doi.org/10.1016/S0167-7152(01)00124-9
  * Zahid et al. (2016) Zahid, H. J., Geller, M. J., Fabricant, D. G., & Hwang, H. S. 2016, ApJ, 832, 203, doi: 10.3847/0004-637X/832/2/203
  * Zahid et al. (2018) Zahid, H. J., Sohn, J., & Geller, M. J. 2018, ApJ, 859, 96, doi: 10.3847/1538-4357/aabe31
  * Zhang et al. (2019) Zhang, H., Primack, J. R., Faber, S. M., et al. 2019, MNRAS, 484, 5170, doi: 10.1093/mnras/stz339
  * Zhuang & Shen (2024) Zhuang, M.-Y., & Shen, Y. 2024, ApJ, 962, 139, doi: 10.3847/1538-4357/ad1183
  * Zonoozi et al. (2019) Zonoozi, A. H., Mahani, H., & Kroupa, P. 2019, MNRAS, 483, 46, doi: 10.1093/mnras/sty2812

##  Appendix A Investigation of Possible PSF Systematics

In this Appendix we show that PSF systematics are unlikely to explain our
results, which are at least an order of magnitude above the conventional weak
lensing regime. Our galaxy shape measurements were done with Galfit (Peng et
al., 2002, McGrath et al., in prep.) assuming a single-component Sérsic model
convolved with an empirical NIRCam PSF. The PSFs are “global” averages based
on stacking all stars that fall within the CEERS footprint in each filter as
described in section 3.2 of Finkelstein et al. (2023). It is customary to
summarize a PSF by only its FWHM, enclosed flux within some aperture and
corresponding limiting magnitude for a point source, but for our purposes, we
also need to estimate its complex ellipticity. Furthermore, in addition to the
global stacked PSF, we also need to consider spatial variations by fitting
individual stars throughout the survey. We use a similar list of stars in each
filter throughout the CEERS footprint as Finkelstein et al. (2023). These
stars were selected via custom cuts on half-light radius and magnitude. For
each star in each filter, we made a $101\times 101$ pixel cutout and discarded
objects with obvious issues (neighbor, detector edge, etc.). This leaves us
with $111-134$ stars per filter, i.e., $\sim 1$ star arcmin-2.

Following standard practice in weak lensing studies (e.g., Kaiser et al.,
1995; Schneider, 2005; Heymans et al., 2005, 2006; Massey et al., 2007a; Jee
et al., 2007; Bridle et al., 2010; Mandelbaum, 2018; Finner et al., 2023a, b;
Cha et al., 2024), we measure the “quadrupole moments” of every individual
star as well as the global stacked PSF image from Finkelstein et al. (2023).
These quadrupole moments characterize the variance and covariance of the light
intensity $I(x,y)$ along the $x$ and $y$ directions. To suppress noise and
ensure that only the “core” of the PSF contributes to the measurements, we
adopt the common choice of a circular Gaussian with $\sigma=3$ pixels as a
weight function $\mathcal{W}(x,y)$ (our results are insensitive to reasonable
variations in $\sigma$). This lets us compute the first moment of the 2D light
distribution which is the centroid ($\bar{x}$,$\bar{y}$):

$\bar{x}=\frac{\sum x\mathcal{W}(x,y)I(x,y)}{\sum\mathcal{W}(x,y)I(x,y)}\;,$
(A1)

$\bar{y}=\frac{\sum y\mathcal{W}(x,y)I(x,y)}{\sum\mathcal{W}(x,y)I(x,y)}\;.$
(A2)

Then we compute the second moments using the pixel indices shifted relative to
the centroid:

$Q_{xx}=\frac{\sum(x-\bar{x})^{2}\mathcal{W}(x,y)I(x,y)}{\sum\mathcal{W}(x,y)I(x,y)}\;,$
(A3)

$Q_{yy}=\frac{\sum(y-\bar{y})^{2}\mathcal{W}(x,y)I(x,y)}{\sum\mathcal{W}(x,y)I(x,y)}\;,$
(A4)

$Q_{xy}=\frac{\sum(x-\bar{x})(y-\bar{y})\mathcal{W}(x,y)I(x,y)}{\sum\mathcal{W}(x,y)I(x,y)}\;.$
(A5)

Finally, this lets us assign to every star (i.e., PSF) a complex ellipticity

$e_{1}+ie_{2}=\frac{Q_{xx}-Q_{yy}+2iQ_{xy}}{Q_{xx}+Q_{yy}+2\sqrt{Q_{xx}Q_{yy}-Q_{xy}^{2}}}\;.$
(A6)

This is analogous to the complex ellipticity that we assigned to every galaxy
in Equation 4 except here we are using quadrupole moments. The complex
magnitude $e=\sqrt{e_{1}^{2}+e_{2}^{2}}$ gives us the PSF ellipticity and the
complex phase $\phi=0.5\rm{arctan2}(e_{2},e_{1})$ gives us the orientation.
The size is defined as $R=\sqrt{Q_{xx}+Q_{yy}}$.

As an example, Figure 19 shows the global empirical PSF in each filter from
Finkelstein et al. (2023) along with our quadrupole moments. As expected, the
ellipticities are very small with $e\lesssim 0.02$. The size $R$ translated to
a Gaussian FWHM, as is common practice, yields similar values as Finkelstein
et al. (2023) with the resolution decreasing from $\sim 0.08"$ in F115W to
$\sim 0.13"$ in F444W. Results for individual stars are similar. Before
proceeding, we caution that although this kind of quadrupole analysis is
standard in weak lensing, the NIRCam PSF is clearly very complicated (i.e.,
not a simple Gaussian or Airy disk). The features outside the “core” of the
PSF are due to JWST’s hexagonally segmented primary mirror and the “tripod”
structure that supports the secondary mirror (Rigby et al., 2023). Other
recent JWST studies (Finner et al., 2023a; Cha et al., 2024; Zhuang & Shen,
2024) also use quadrupole moments or fit an elliptical 2D Gaussian to the
“core” of the PSF and arrive at similar conclusions as us, but we recommend
more detailed follow-up analysis.

Figure 19: The global empirical PSF from stacking stars in each filter from
Finkelstein et al. (2023) along with our quadrupole moment measurements. The
white ellipse illustrates the ellipticity, orientation and size of the
quadrupole moments that characterize the core. Naturally, the ellipticities
are small with $e\lesssim 0.02$ so these ellipses are nearly circular. All
images are on a log scale to emphasize JWST’s unique diffraction features
outside the “core” of the PSF. On a linear scale, the size $R$ derived from
quadrupole moments adequately reflects the extent of the core. Results for
individual stars are similar.

Figure 20 shows the distribution of our individual stars and the global PSF in
the complex $(e_{1},e_{2})$ plane. Nearly all stars in every filter have
$e<0.05$, i.e., the PSF is very round as expected. The short-wavelength NIRCam
filters (F115W, F150W, F200W) have average components $\langle
e_{1}\rangle\sim 0.02$ and $\langle e_{2}\rangle\sim-0.01$ and there is
clearly a bias towards the lower-right quadrant. The long-wavelength NIRCam
filters (F277W, F356W, F444W) have substantially rounder PSFs as evidenced by
their smaller $\langle e_{1}\rangle\lesssim 0.01$ and $\langle
e_{2}\rangle\sim\pm 0.001$. These percent-level PSF systematics are unlikely
to explain our much larger observed $\langle e\rangle\gtrsim 0.1$ when
averaging over galaxies in individual NIRCam chips, modules and pointings
(Figure 18). At most, they would lead to an additive bias in $\langle
e\rangle$ of a few percent, not tens of percent.

Figure 20: Distribution of our individual stars (blue points) and global
stacked PSF (yellow star) in the complex $(e_{1},e_{2})$ plane. This is in the
standard “north up” frame. Nearly all stars in every filter have $e<0.05$,
i.e., the PSF is very round. It is slightly more elliptical in the short-
wavelength filters (F115W, F150W, F200W) with $\langle e_{1}\rangle\sim 0.02$
and $\langle e_{2}\rangle\sim-0.01$ and there is clearly a bias towards the
lower-right quadrant. The long-wavelength filters (F277W, F356W, F444W) have
even rounder PSFs as expected. These percent-level PSF systematics cannot
explain our much larger $\langle e\rangle\gtrsim 10\%$ when averaging over
background galaxies (Figure 18) but must be corrected in any future analysis
of the conventional weak lensing regime.

Figure 21 shows a map of the PSF ellipticity and orientation across the entire
CEERS footprint using our individual stars. Most stars have very small
ellipticity but there are a handful with $e\sim 0.05-0.08$. The cores of these
stars do look slightly elliptical and it is not clear whether they are
outliers since other stars near them show smaller ellipticities. The
orientations of PSF stars in the long-wavelength filters appear random due to
their very small ellipticities. In contrast, for the short-wavelength filters,
there is an overall pattern of many stars pointing towards the upper-left of
each panel ($45^{\circ}$ east of north). However, again, since the individual
PSF ellipticities are very small, their additive bias due to these alignments
cannot contribute more than a few percent to our much larger observed $\langle
e\rangle>10\%$ from averaging over background galaxies.

Figure 21: Ellipticity and orientation of the PSF at the positions of
individual stars across the entire CEERS footprint. For convenience, the
ellipticity is indicated both with the line length (scale bar in lower left)
and with the colorbar. The ellipticity and orientation of the stacked PSF is
shown in the left center of each panel. This is in the standard “north up”
frame. The ellipticities are very small and appear random in the long-
wavelength filters. In contrast, for the short-wavelength filters, there are
systematic PSF alignments towards $45^{\circ}$ east of north but the small
ellipticities mean that these will cause an additive bias of at most a few
percent which is an order of magnitude lower than our $\langle e\rangle\gtrsim
10\%$ from averaging background galaxies (Figure 18). However, future studies
in the conventional weak lensing regime must correct for this bias.

Finally, Figure 22 shows that our actual background galaxy sample spans a much
larger range in this same complex $(e_{1},e_{2})$ plane. Only $\sim 2\%$ of
our galaxies have $e<0.05$ similar to our PSF stars and they are roughly
uniformly split between the different quadrants (phase angles). Note that the
majority of our galaxies ($2972/3848\sim 77\%$) are assigned to a short-
wavelength NIRCam filter (F115W, F150W, F200W) to obtain the rest-frame
morphology corresponding to their redshift. Moreover, our smallest size
galaxies (which would be preferentially subject to PSF systematics) are also
roughly uniformly spread out in this plane, i.e., not biased towards any one
quadrant. This is true regardless of whether we look at the combined sample or
only galaxies in individual pointings. Given that the PSF tends to round out
objects, the fact that our galaxies have intrinsic ellipticities that are much
larger than that of the stars implies that correcting for the small bias
towards the lower-right quadrant in Figure 20 would not cause most of our
galaxies to change quadrants. Thus we do not believe that the percent-level
PSF systematics remaining in the data can explain our observed alignments in
regions with $\langle e\rangle\gtrsim 10\%$. However, measuring the precise
amplitude of the shear variance $\langle\overline{\gamma}^{2}\rangle$ will
require correcting galaxy shapes for this spatially-varying PSF because the
expected signal is weak, especially on larger scales (Figure 18). It is
unclear whether CEERS itself has enough stars to interpolate the PSF shape to
the position of every galaxy so a more sophisticated strategy may be required.
We defer that to a future more detailed weak lensing analysis of “blank” JWST
deep fields.

Figure 22: Distribution of our background galaxy sample in the complex
$(e_{1},e_{2})$ plane (standard “north-up” frame). The large center panel
shows our combined background sample whereas the smaller panels show the
distributions in individual pointings. In all cases, $\sim 2\%$ of our
galaxies have $e<0.05$ which is the PSF-dominated regime (within the central
black circle). Our galaxies also tend to be roughly uniformly distributed
between the four quadrants. This is particularly true for our smaller-size
galaxies (bluer colored points) which would be most susceptible to PSF
effects. The fact that most of our galaxies have intrinsic ellipticities that
are much larger than that of our stars implies that the small bias towards the
lower-right quadrant in Figure 20 is unlikely to dictate the orientations and
hence alignments of our galaxies.