# Prospects for studying the mass and gas in protoclusters with future CMB
observations

Anna Gardner,1 Eric Baxter,1 Srinivasan Raghunathan,2 Weiguang Cui,3 Daniel
Ceverino4,5
1Institute for Astronomy, University of Hawai‘i, 2680 Woodlawn Drive,
Honolulu, HI 96822, USA
2Center for AstroPhysical Surveys, National Center for Supercomputing
Applications, Urbana, IL 61801, USA
3Institute for Astronomy, University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, United Kingdom
4Universidad Autonoma de Madrid, Ciudad Universitaria de Cantoblanco, E-28049
Madrid, Spain
5CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, E-28049 Madrid,
Spain

###### Abstract

Protoclusters are the progenitors of massive galaxy clusters. Understanding
the properties of these structures is important for building a complete
picture of cluster formation and for understanding the impact of environment
on galaxy evolution. Future cosmic microwave background (CMB) surveys may
provide insight into the properties of protoclusters via observations of the
thermal Sunyaev Zel’dovich (SZ) effect and gravitational lensing. Using
realistic hydrodynamical simulations of protoclusters from the The Three
Hundred Project, we forecast the ability of CMB Stage 4-like (CMB-S4)
experiments to detect and characterize protoclusters with observations of
these two signals. For protoclusters that are the progenitors of clusters at
$z=0$ with $M_{200c}\gtrsim 10^{15}\,M_{\odot}$ we find that the S4-Ultra deep
survey has a roughly 20% chance of detecting the main halos in these
structures with ${\rm SNR}>5$ at $z\sim 2$ and a 10% chance of detecting them
at $z\sim 2.5$, where these probabilities include the impacts of noise, CMB
foregrounds, and the different possible evolutionary histories of the
structures. On the other hand, if protoclusters can be identified using
alternative means, such as via galaxy surveys like LSST and Euclid, CMB-S4
will be able to obtain high signal-to-noise measurements of their stacked
lensing and SZ signals, providing a way to measure their average mass and gas
content. With a sample of 2700 protoclusters at $z=3$, the CMB-S4 wide survey
can measure the stacked SZ signal with a signal-to-noise of 7.2, and the
stacked lensing signal with a signal-to-noise of 5.7. Future CMB surveys thus
offer exciting prospects for understanding the properties of protoclusters.

###### keywords:

galaxies: clusters: general – large-scale structure of Universe – galaxies:
clusters: intracluster medium – cosmic background radiation

††pagerange: Prospects for studying the mass and gas in protoclusters with
future CMB observations–Acknowledgments

##  1 Introduction

Galaxy clusters are the most massive gravitationally bound structures in the
Universe, with masses of roughly $10^{14}$ to $10^{15}\,M_{\odot}$ and virial
radii of order 1 Mpc. Although the definition of a protocluster is not
entirely consistent across the literature, one common and practical definition
is that a protocluster is a structure that will collapse into a galaxy cluster
by $z\geq 0$ (see also Muldrew:2015). By studying protoclusters, we can learn
about the evolutionary histories of galaxy clusters, explore the impact of
environment on galaxy formation (e.g. Steidel:2005; Muldrew:2017), learn about
processes driving reionization (e.g. Chiang_2017), and more (for a review see
Overzier:2016).

Protoclusters are difficult to detect because they are at high redshift,
[typically reside in halos that are less massive than virialized clusters at
low redshift, and they can extend over tens of comoving Mpc. Protoclusters are
also rare, having formed from extreme density fluctuations. Detecting
protoclusters therefore requires wide-field surveys in addition to high
sensitivity.] For these reasons, the number of protoclusters that have been
detected to date is small, on the order of 100 (and of order 10 above $z\sim
4$). Several techniques have been used to detect protoclusters, including
looking for overdensities of galaxies selected via e.g. dropout methods in
photometric surveys (e.g. Toshikawa:2018), looking for overdensities of star-
forming galaxies at high redshift via their redshifted emission
(Clements:2014; miller18, for example), and looking for correlated
Lyman-$\alpha$ absorption along lines of sight to background galaxies
(Newman_2022).

The majority of protoclusters found to date have been identified via galaxy
overdensities. However, protocluster selection via galaxy overdensities
depends strongly on the galaxy sample, and the identified overdensities may or
may not collapse into massive clusters at $z=0$ (Cui2020). [Moreover, the
completeness and purity of protocluster samples identified via galaxy
overdensities depend strongly on e.g. the choice of aperture used to identify
the overdensities (Muldrew:2015; Lovell:2018). The purity and completeness of
protocluster regions identified via galaxy density has also been investigated
using The300 simulations in Cui2020]. The galaxy populations within
protoclusters are often used to form rough mass estimates of these structures
by, for instance, using stellar mass–halo mass relations (e.g. Laporte:2022)
or by relating the galaxy overdensity to the underlying mass overdensity by
assuming a value for the galaxy bias (e.g Toshikawa:2012). These estimates
typically carry large uncertainties and often rely on assumptions about the
member galaxy properties.

The aim of this work is to investigate the prospects for detecting and
characterizing the properties of protoclusters without using their member
galaxies. Instead, we consider how the cosmic microwave background (CMB) can
be used to study protoclusters. Because the CMB originates from very high
redshift ($z\sim 1089$), it provides a backlight to all clusters and
protoclusters in the Universe. CMB photons travelling near to clusters and
protoclusters can be perturbed by these structures in several different ways.
We focus on two such effects here. First, CMB photons can inverse Compton
scatter with populations of thermal electrons in protoclusters, leading to the
thermal Sunyaev Zel’dovich (SZ) effect (Sunyaev:1972). The amplitude of this
effect is proportional to the Compton-$y$ parameter, which is sensitive to the
line-of-sight integral of the electron gas pressure. The SZ effect is
routinely used to detect galaxy clusters with current CMB data, providing
catalogs of clusters with masses $M\gtrsim{\rm few}\times 10^{14}\,M_{\odot}$
out to redshifts of $z\sim 1.5$ (e.g. Bleem:2015; Hilton:2021;
Planck:clusters). Since protoclusters are expected to have significant ionized
gas content, they may also produce detectable SZ signatures. However,
detecting the SZ signals from protoclusters will be challenging for several
reasons. For one, the amount of mass in a protocluster that is contained in a
single virialized halo is significantly less than the mass of the low-redshift
cluster into which the structure will evolve (Muldrew:2015). As a result, the
SZ signals from protoclusters are expected to be much lower than those of low-
redshift clusters. The expected scaling of the SZ signal, $Y$, with halo mass,
$M$, is $Y\sim M^{5/3}$ for self-similar halos (Kaiser:1986) [(see e.g.
Pratt:2020 for recent observational confirmation of this prediction)]. Since
the typical cluster with $M\sim 10^{15}\,M_{\odot}$ at $z=0$ grew from a
protocluster with a main halo of mass $\sim 10^{13.6}$ at $z=2$ (see Fig. 1),
the average protocluster at $z=2$ will have an SZ signal which is roughly a
factor of 200 lower than that of a typical $z=0$ cluster. [Since typical wide-
field SZ surveys like those from the South Pole Telescope (SPT; Bleem:2015),
the Atacama Cosmology Telescope (ACT; Hilton:2021) and Planck
(Planckszclusters) detect clusters with signal-to-noise of order tens, it is
not too surprising that detections of many protoclusters have not been
reported by these surveys. ] This calculation ignores mass outside of the main
halo, but doing so is reasonable on account of the $M^{5/3}$ scaling of the SZ
signal. Since protocluster halos may not be self-similar, this scaling is also
not expected to hold exactly, but should give a reasonable indication of the
challenge associated with high-redshift protocluster detection via the SZ (see
Sembolini2014, for a prediction of steeper slope). Moreover, we expect halos
at high redshift to be less concentrated than their low-redshift counterparts,
leading to additional suppression of the SZ signal.111Halo concentration is
positively correlated with SZ signal at fixed halo mass, see e.g. Baxter:2023.
Despite these challenges, detecting the SZ signals from protoclusters is an
exciting prospect, as such a detection would provide a window into the gas
content and thermal properties of these structures. Indeed, during the
completion of this work, the detection of the SZ signal from a protocluster
was reported by Mascolo:2023.

A second way that clusters and protoclusters may leave an imprint on the CMB
is via gravitational lensing. The paths of CMB photons travelling near massive
objects like clusters or protoclusters will be perturbed, leading to an
observable impact on the CMB (Seljak:2000; Dodelson:2004; Lewis:2006; hu07).
Gravitational lensing of the CMB by galaxy clusters can now be detected and
used to calibrate cluster mass-observable relations (madhavacheril15;
Baxter:2015; Baxter:2018; raghunathan19; raghunathan19c; madhavacheril20). We
consider here the possibility of using CMB lensing to constrain the masses of
protoclusters. Gravitational lensing in general is a powerful tool for
measuring the masses of (proto)clusters because it is sensitive to all forms
of mass, including dark matter. Mass estimation methods that rely on visible
matter, on the other hand, must make assumptions about the connection between
the visible matter and the underlying dark matter, which makes up most of the
mass in these structures (e.g. that the visible and dark matter are in virial
equilibrium). Such assumptions are likely to be particularly bad for
protoclusters, which may still be in the process of forming, and which may
contain galaxies that differ significantly from those at low redshifts or in
the field. CMB lensing has the potential to be a particularly powerful tool
for measuring protocluster masses because — unlike galaxy lensing — it can be
measured around structures at very high redshift. Gravitational lensing
measurements relying on galaxy lensing, on the other hand, will be challenging
to apply to protoclusters because they will require measurements of the shapes
of background galaxies at higher redshifts than the protoclusters
(Madhavacheril:2017).

In this work, we use realistic hydrodynamical simulations of protoclusters
from The Three Hundred Project222http://the300-project.org/ (The300 for short;
Cui2018) to investigate the extent to which future CMB surveys can measure the
SZ and CMB lensing signals from protoclusters. Our use of hydrodynamical
simulations to make these forecasts is novel and important because actively
forming protoclusters may not be described by simple analytical models.
Previous work making similar forecasts (e.g raghunathan22) has generally
relied on analytical models developed for virialized clusters in order to make
forecasts for protocluster measurements. [Moreover, with these simulations, we
can explicitly investigate the impact of each protocluster’s evolutionary
history on its SZ and lensing signals. We will find that there is significant
scatter in these signals as a result of different evolutionary histories, an
effect which is often missed in analytic forecasts for protocluster
measurements.] While some details of protocluster evolution (such as the
impact of AGN feedback and environment-dependent quenching of star formation)
remain poorly understood and thus unlikely to be fully captured in the
simulations, we expect that protocluster mass distributions and bulk gas
properties will be accurately modeled. It is these quantities which are
relevant to determining the SZ and CMB lensing signals. In any case, our
intent here is to make a forecast for the expected observational signatures of
protoclusters given a specific model (i.e. the Gadget-X model in The300
simulations); departures from this model in real protoclusters would be
interesting.

Near-term CMB surveys including Simons Observatory (SimonsObservatory) and CMB
Stage 4 (CMBS4) will map large areas of the microwave sky to unprecedented
depth. We will show that these observations can place interesting constraints
on the SZ and CMB lensing signatures of protoclusters. We consider two
possibilities: (1) detecting protoclusters with the CMB survey and
characterizing their individual properties, and (2) a “stacking” approach,
where some other survey — such as a galaxy survey — is used to identify
protocluster locations, and the CMB survey is then used to constrain the
average mass and gas properties of the resultant samples.

As mentioned above, CMB surveys like Planck have already been used to detect
protoclusters by searching for emission from protocluster member galaxies
(Clements:2014). These detections rely on the high-frequency channels of
Planck (353 to 857 GHz) to detect infrared emission from dusty, star-forming
galaxies that has been highly redshifted. Here, in contrast, we focus on
detecting signals from the impact of protoclusters on CMB photons, rather than
emission directly from the protoclusters themselves. These “indirect” signals
have the advantage that they are largely independent of assumptions about the
protocluster member galaxies.

The organization of the paper is as follows: In §2 we discuss the simulation
data products that are used in this work and our analysis of these data
products, in §3 we present our main results, and we conclude in §4. Throughout
this work we adopt a flat $\Lambda$ and cold dark matter cosmological model
with parameters consistent with those of The300 (Cui2018), which are in turn
consistent with the results of the Planck mission: $h=0.678$, $\Omega_{\rm
m}=0.307$, $\Omega_{\rm b}=0.048$, $n_{s}=0.96$, $\sigma_{8}=0.8228$
(planck20_cosmo). Unless specified otherwise, halo mass refers to $M_{200c}$,
i.e. the mass contained within a sphere of radius $R_{200c}$ centered on the
halo, where $R_{200c}$ is the radius at which the mean enclosed density is 200
times the critical density of the Universe at the redshift of the halo.

##  2 Simulated protocluster observations

Figure 1: The redshift evolution of the mass distribution of the main
progenitor halos of clusters in The300 simulations, which are selected based
on a mass cut of $M\gtrsim 10^{15}M_{\odot}$ at $z=0$. At high redshift, the
main halo in a protocluster may constitute only a small fraction of the
structure’s total mass. [The numbers in parentheses in the legend indicate the
mean and standard deviation of $M_{200c}$ for each redshift, in units of
$10^{15}\,M_{\odot}$]

###  2.1 The300 project

The300 Project (Cui2020) is a re-simulation of a mass-complete sample of 324
galaxy clusters from the MultiDark Planck 2 (MDPL2) $N$-body simulation
(Klypin2016)333https://www.cosmosim.org/cms/simulations/mdpl2 with $1{\rm
Gpc}/h$ simulation box size (see Zhang2022; deAndres2022NatAs; deAndres2023,
for the benefit of this particular setup). These clusters are identified with
the rockstar (Rockstar) halo finder. A large high-resolution region $\sim
15{\rm Mpc}/h$ in radius (see Rost2021; Kuchner2020; Kuchner2021; Kuchner2022,
for the studies of filaments around clusters thanks to this large region) at
$z=0$ surrounding the central cluster is traced back to a very high redshift
$z=120$ and the initial conditions (ICs) are generated by the parallel
GINNUNGAGAP code444https://github.com/ginnungagapgroup/ginnungagap including
gas particles in the high-resolution region by splitting the dark matter
particles and using low-resolution particles to sample the tidal fields
outside of the high-resolution region in multiple layers. From these ICs,
three hydrodynamic simulation codes are used to evolve the clusters: Gadget-
MUSIC (Sembolini2013); Gadget-X (Rasia2015); Gizmo-SIMBA (Dave2016; Dave2019;
Cui2022), [all of which are updated versions of Gadget2 (Springel2005), a
widely used N-body/SPH simulation code]. [The three codes include models for
baryonic processes such as gas cooling to form stars and the enrichment of gas
by supernovae. Gadget-X and Gizmo-SIMBAalso include black hole models and
prescriptions for feedback from active galactic nuclei. We refer readers
Cui2022 for more detailed comparisons.]

In this work, we use only the simulations from the Gadget-X run, which agrees
best with observed gas properties (Li2021; Li2023).

The mass distribution of the main halos in The300 as a function of redshift is
shown in Fig. 1. At $z=0$, the distribution corresponds to a mass cut at
$M\gtrsim 10^{15}M_{\odot}$. At higher redshifts, the mass distribution of the
main halos becomes broader and shifts to lower values. Note that at high
redshift, the main halo in a protocluster may constitute only a small fraction
of the total mass of the structure (Muldrew:2015). Furthermore, it is worth
noting that although the sample is selected to be mass-complete at $z=0$,
halos in The300 are not always the most massive halos at high redshift (see
Cui2018; Cui2022EPJWC, for the mass completeness at different redshifts in
these simulations).

###  2.2 SZ and lensing maps

To characterize the SZ signals from protoclusters we produce simulated Compton
$y$ maps at several snapshots from The300 simulations. Compton $y$ is related
to the line-of-sight integral of the electron gas pressure, $P_{e}$, via

$y=\frac{\sigma_{T}}{m_{e}c^{2}}\int dl\,P_{e}(l),$ (1)

where $l$ is the line-of-sight distance. We use $Y$ to represent the integral
of $y$ over some circular aperture centered on the cluster:

$Y=2\pi\int_{0}^{R_{\rm max}}r\,dr\,y(r),$ (2)

where $r$ is the projected distance from the cluster center and $R_{\rm max}$
is the size of the aperture. Following common convention, we set $R_{\rm
max}=R_{200c}$. The observed, frequency-dependent intensity change in the CMB
maps due to the SZ effect is then proportional to Compton $y$
(Carlstrom:2002). [We will work entirely in Compton $y$ in this analysis,
rather than frequency-dependent intensity. In principle, there may be some
advantages to working with the raw intensity maps, particularly with regard to
suppressing contamination from foregrounds (e.g. Chiang:2020). However, our
choice greatly simplifies the analysis, and we expect it to have little impact
on the forecast signal to noise, which is the main focus of this work. For the
protoclusters that we focus on, their masses are sufficiently small that the
relativistic thermal SZ effect is expected to be negligible. The kinematic SZ
effect can also be ignored given its small amplitude relative to the thermal
SZ.]

In detail, these maps, which cover almost the full high-resolution region of
the simulation (24 comoving ${\rm Mpc}/h$ on a side), are generated by the
PyMSZ package555https://github.com/weiguangcui/pymsz (see Cui2018 for a more
in depth description of the process of generating the mock $y$-maps). The
angular resolution of the mock maps is initially fixed to 3$"$, i.e. much
higher resolution than expected for the future CMB surveys of interest; we
discuss the degradation of the resolution due to the telescope beam in §2.3.
Because the angular scale is fixed, the mock maps have different numbers of
pixels and pixel sizes to cover the 24 ${\rm Mpc}/h$ on a side region at the
different redshifts. The line-of-sight depth used when computing $y$ is also
24 comoving ${\rm Mpc}/h$; we discuss our inclusion of noise from foreground
structure at larger distances in §2.3. The mock maps are always centred on the
centre of the whole simulation box, not at the most massive halos, which can
be very far away from the simulation centre at high redshift. Though at low
redshift, this region is much larger than the halo’s virial radius, this
guarantees that the whole protocluster region at high redshift is included.

Figure 2: Power spectra of the residual noise and foregrounds after applying
the minimum variance ILC weights to data from the S4-Wide (left panel) and
S4-Ultra deep (right panel) surveys are shown using the thick black curves.
The thin coloured curves correspond to experimental noise in individual
frequency bands. The lensed CMB power spectrum $C_{\ell}^{TT}$ is shown as the
gray solid curves for reference. Figure 3: The total noise and foreground
residuals expected in the MV Compton-$y$ map for S4-Wide (green) and S4-Ultra
deep (red) surveys. The gray line indicates the Compton-$y$ power spectrum
from large scale structure, which acts as an additional source of noise for
the protocluster measurements.  Figure 4: CMB lensing noise power spectra
considered in this work. The power spectra here are given in terms of the
lensing potential $\phi$ rather than the convergence $\kappa$. The green (red)
curves represent the noise and foreground contributions for the S4-Wide
(S4-Ultra deep) surveys, while the grey curve represents the lensing power
spectrum from late-time large scale structure, which constitutes a source of
noise for measurement of the protocluster lensing signal.

The lensing maps are produced in a similar fashion to the SZ maps. Under the
Born approximation, the lensing deflection can be computed by integrating
along the paths of undeflected photons. In this case, the lensing convergence
is related to the surface density in the plane of the lens, $\Sigma$, via

$\kappa=\frac{\Sigma}{\Sigma_{\rm crit}},$ (3)

where $\Sigma_{\rm crit}$ is the critical surface density given by

$\Sigma_{\rm crit}=\frac{c^{2}D_{s}}{4\pi GD_{l}D_{ls}},$ (4)

where $D_{s}$ is the angular diameter distance to the source (in this case the
last scattering surface, where the CMB originates), $D_{l}$ is the angular
diameter distance to the lens, and $D_{ls}$ is the angular diameter distance
between them. We compute the surface density by integrating the simulation
particle data along a randomly chosen line of sight.

Note that $\kappa$ is not directly observable; rather, one measures the CMB
temperature and polarization signals across the sky, and uses the statistical
properties of these fluctuations to estimate $\kappa$ (e.g. hu02). Here, for
simplicity we will work entirely with $\kappa$ rather than temperature or
polarization maps.

The SZ and lensing maps are generated at several redshifts:
$z\in[3.5,3,2.5,2.0,0.9,0.02]$ spanning the most interesting regime for
detecting and characterizing protoclusters.

###  2.3 Mock observations

In order to make forecasts for the ability of future CMB surveys to constrain
the properties of protoclusters, we simulate the impact of noise,
astrophysical foregrounds, and the telescope beam on the SZ and lensing
observations. We consider both the wide (S4-Wide) and deep (S4-Ultra deep)
surveys from the CMB-S4 experiment (cmbs4-sb1; cmbs4collab19).

In addition to noise, CMB foregrounds such as the cosmic infrared background
and radio sources present a significant challenge for estimation of Compton
$y$ and CMB temperature/polarization (i.e. the relevant quantities for
measuring CMB lensing). To reduce the impact of foregrounds, one must rely on
some form of component separation. [For our forecasts, we consider the use of
internal linear combination (ILC) methods to estimate $y$ and the CMB
temperature signal (e.g Tegmark:2003; Eriksen:2004). Such methods form a
linear combination of the observations at different frequencies to estimate
the component of choice. Here, we choose weights in the linear combination
such that there is unit response to the component of interest, i.e. CMB
temperature (for the lensing forecasts) or $y$ (for the SZ forecasts) and the
variance is minimized. Note that for the SZ forecasts, the primary CMB
fluctuations constitute a noise source. We refer to this set of ILC weights as
minimum variance or (MV).] This procedure accounts for additional variance
introduced by foregrounds, but does not account for possible biases due to
foregrounds, as we discuss below.

We model the foregrounds using measurements made by SPT reichardt21). We
include signals from primary CMB, kinematic SZ, radio galaxies and dusty star
forming galaxies. The noise is decomposed into atmospheric and instrumental
parts. The resultant temperature power spectra are shown in Fig. 2. For more
details about our noise and foreground modeling, we refer readers to
raghunathan22; raghunathan22b. In Fig. 3, we show the combined noise and
foreground power spectra in the Compton-$y$ map for S4-Wide in green and
S4-Ultra deep in red.

Since Compton-$y$ is sensitive to the integrated pressure along the line of
sight, large scale structure (LSS) will also contribute to the total variance
when studying the signals from protoclusters. The grey curve in the figure
represents the calculation of the power spectrum of Compton $y$ from LSS
including halos greater than $M^{10}\,M_{\odot}$.666Our calculation of the LSS
contribution includes only the one-halo term, which is expected to dominate
across the scales of interest (Komatsu:2002). The total noise power used in
our analysis is given by the sum of the LSS and residual noise and foreground
contributions. In principle, the LSS contribution to the noise could be
reduced by identifying and masking halos with large SZ signals that are
unassociated with protoclusters; by including the contributions from all halos
above a low threshold we are therefore being conservative. Note that the noise
spectra shown in Fig. 3 and Fig. 4 do not include the impact of the
instrumental beam, which will suppress the noise at high $\ell$; we discuss
the inclusion of the beam below.

[Our estimate of the CMB lensing noise assumes that quadratic estimator (QE)
methods (hu02; okamoto03; hu07) are used to estimate the CMB lensing signal.
The QE extracts the lensing signal using the lensing-induced correlations
between different multipoles of the temperature and polarization fields (which
are expected to be zero for a Gaussian and isotropic primordial CMB in the
absence of lensing). Chance correlations between different multipoles of the
fluctuating CMB field, as well as instrumental noise, then contribute to the
total CMB lensing noise. We estimate the noise power of the CMB lensing maps
following hu02.]

The presence of non-Gaussian foregrounds can also introduce such mode
coupling, leading to biases in the reconstructed lensing signal. Particularly
worrying is the SZ signal from the (proto)clusters themselves, which can bias
the lensing estimates for low-redshift clusters (Baxter:2015; Baxter:2018).
However, this bias is likely small for protoclusters given the expected
$M^{5/3}$ scaling of the SZ signal. Additional potential sources of bias
include emission from galaxies residing within the protoclusters. In general,
recently developed methods may be able to mitigate some of these foreground
biases with negligible impact to the signal-to-noise (madhavacheril18;
raghunathan19b; raghunathan19c; Hadzhiyska:2019; Horowitz:2019; levy23). Our
CMB lensing noise estimate is derived assuming a combination of temperature
and polarization measurements, and since the foreground signals are largely
unpolarized, the lensing measurements based on polarization data are expected
to be immune to foreground biases. For these reasons, we will ignore potential
biases in the lensing signal due to foregrounds, postponing a detailed study
to future work. [Our final estimate of the CMB lensing noise power spectrum is
shown in Fig. 4.]

[Analogously to the SZ forecasts, we must also include the impact of CMB
lensing fluctuations due to large scale structure, which constitute a noise
source for our lensing forecasts. We compute this contribution using
CAMB777camb.info (Lewis:1999bs; Howlett:2012mh). This contribution is shown as
the grey curve in Fig. 4.]

[Throughout this analysis, we assume that all sources of noise for our
forecasts are statistically isotropic on the sky. This is an excellent
approximation for the extragalactic foregrounds. Moreover, across each survey
region (S4-Wide and S4-Ultra deep), CMB-S4 instrumental noise is expected to
be close to uniform. Galactic foregrounds, on the other hand, may introduce
anisotropy into the noise. We note, though, that our forecasts assume
conservative values of sky coverage to account for removing regions of high
galactic emission near the galactic plane. Outside of these regions, the
contributions of galactic foregrounds to the total foreground power are
subdominant. ]

We generate Gaussian realizations of the SZ and lensing noise power spectra at
the same pixel scale as the simulated SZ and lensing maps. While the true
noise is unlikely to be purely Gaussian,888For instance, the “noise”
contributed by LSS must be non-Gaussian at some level. this simple approach
should be sufficient for the purposes of our forecasts. To account for the
impact of the telescope beam on the measurements, we apply Gaussian smoothing
to the signal and noise maps999In reality, instrumental noise is not impacted
by the telescope beam. However, in a real-space analysis like ours, one would
typically smooth the noisy CMB lensing maps with an effective beam to suppress
small-scale noise. Our application of the beam to the noise maps reflects this
procedure.  with $\theta_{\rm FWHM}=1^{\prime}$, close to the expected
effective resolution of the CMB lensing maps from CMB-S4. The smoothed signal
and noise maps are then summed to obtain the final mock-observed protocluster
maps. [For essentially all of the protoclusters that we consider, the assumed
beam size is smaller than the protocluster virial radii. Consequently, we do
not expect higher resolution observations to yield much improvement in signal-
to-noise per protocluster (although higher resolution observations might be
useful for, e.g., distinguishing features in protocluster SZ profiles). ]

###  2.4 SZ and lensing profile measurements

A powerful way to characterize the lensing and SZ signals from protoclusters
will be to measure their radial profiles. We compute such profiles for each of
the 324 simulated clusters by measuring the average of the lensing and SZ
values in annuli centered on the protocluster, where the center is defined as
the location of the maximum in the SZ or lensing maps. At the resolutions
accessible to future wide-field CMB experiments, the distinction between these
two choices is unimportant. The annuli used to measure the profiles are
defined in projected radius. We use five bins from a minimum scale of $R_{\rm
min}=0$ to a maximum scale of $R_{\rm max}=4\,{\rm Mpc}$ for both the SZE and
lensing signals. Note that in a real data analysis, using bins of projected
radius (as opposed to angular bins) would require knowledge of the
protocluster redshift, which could be obtained using e.g. spectroscopic
followup. The CMB lensing and SZ measurements have noise power on large
scales. These large-scale noise fluctuations will contribute to the covariance
of the real-space profile measurements, increasing the diagonals and leading
to large off-diagonals. This is not a problem in principle for an analysis
that includes the full covariance, but it does have the downside of making
plots of the profile measurements much harder to interpret. To reduce this
effect, we subtract off a measurement of the local background from each of the
protocluster profiles. This background estimate is derived from an annulus of
width $0.8\,{\rm Mpc}$, with minimum radius equal to the largest radius
profile measurement.

In order to make forecasts for detectability of the SZ and lensing signals
with CMB-S4, we compute the noise contribution to the profiles using the mock
noise maps described in §2 using the same binning choices and background
annulus subtraction as for the cluster measurements. Since the noise is
uncorrelated with the clusters and is additive, our estimated noise covariance
for the profile of a single cluster is given by the covariance of the noise
profiles.

As noted previously, we also consider a stacking analysis in which we assume
that protocluster positions have been determined via some other survey (e.g. a
galaxy survey) and the signal-to-noise of the SZ and CMB lensing measurements
is then enhanced by averaging across multiple protoclusters (which may or may
not be individually detected in the CMB survey). We form stacked profiles by
averaging the individual protocluster profiles across all of the simulated
protoclusters. Two sources contribute variance to the stacked profile
measurements: noise/foreground fluctuations in the SZ or CMB lensing maps, and
intrinsic variation in the protocluster profiles themselves, both of which are
included in our analysis. The total covariance of the stacked profile
measurements is computed from the sum of the covariances of the noiseless
protocluster profiles and the noise-only covariance. We compute the covariance
of the stacked profile measurements using all clusters in The300, and then re-
scale the errorbars to the expected number of observed clusters (see §2.5).

By chance, there is one very-high mass cluster in our sample that is
essentially fully virialized at high redshift, leading to it having large SZ
and lensing signals (cluster 3 in Fig. 5). This cluster contributes
significantly to the estimated variance of the stacked profile measurements.
We therefore remove this cluster when computing averages across the full
cluster sample. If necessary, a similar approach could also be taken in a real
data analysis, with any anomalously high SNR clusters analyzed separately.

Figure 5: The evolution of three example protoclusters, as seen via the SZ
effect ($y$, blue) and CMB lensing ($\kappa$, red). Each panel shows a region
24 comoving $h^{-1}{\rm Mpc}$ on a side, and different rows correspond to
different redshifts (indicated at left). The color scales are clipped at high
$y$ and $\kappa$.

(a) Sunyaev Zel’dovich effect (b) CMB lensing

Figure 6: The signal-to-noise distribution of the SZ (left) and lensing
(right) detections for individual simulated clusters (colored solid curves).
For reference, the black dashed curve represents the calculated signal-to-
noise distribution when there is no cluster signal (i.e. noise only). Values
greater than zero are expected in this case due to random fluctuations. Beyond
$z\sim 2.5$, some protoclusters will be detected at high significance via
their SZ signatures. On the other hand, essentially no protoclusters can be
detected beyond $z\sim 2$ via their CMB lensing signals.

###  2.5 Expected number of protoclusters

The uncertainties on the stacked profile measurements depend on the number of
protoclusters that we average over. Since we define protoclusters as the
progenitors of a selection of redshift $z=0$ clusters, the number of such
objects per comoving volume remains fixed. The number of protoclusters in the
volume probed by CMB-S4 observations is then given by the number of clusters
in The300 multiplied by the ratio of the comoving volume probed by CMB-S4 to
the volume simulated (1 $({\rm Gpc}/h)^{3}$). We adopt sky coverage of $f_{\rm
sky}=0.5$ for S4-Wide and $f_{\rm sky}=0.03$ for S4-Ultra deep
(cmbs4collab19). We consider redshift bins of width $\Delta z=0.1$ centered on
the redshift values given above. This bin width is small enough that the
evolution of protoclusters across these slices is not dramatic, and large
enough that there is a significant number of protoclusters in each redshift
bin. We find that the number of (proto)clusters in each redshift slice is then
$N_{c}\approx[2600,2700,2600,2500]$ for S4-Wide, and
$N_{c}\approx[160,160,160,150]$ for S4-Ultra deep. Of course, these numbers
will be different for different choices of $\Delta z$; also, assigning
protoclusters to redshift bins will require obtaining redshifts, which may be
obtained via spectroscopic follow up.

Note that this simple approach to estimating the number of protoclusters in
the stack assumes that all protoclusters meeting our definition (i.e. that
they give rise to clusters with $M\gtrsim 10^{15}\,M_{\odot}$ at $z=0$) within
the volume probed by the CMB-S4 observations are detected. This assumption is
likely to be unrealistic in several respects. First, the threshold of
$10^{15}\,M_{\odot}/h$ is largely arbitrary, set by choices made in the design
of The300 simulations. Still, this threshold is a reasonable and conservative
choice in the sense that it corresponds to a very high mass cluster. It is
likely that future galaxy-based protocluster surveys will be able to probe the
progenitors of less massive clusters, given that current surveys can identify
protoclusters expected to be progenitors of clusters with $M\gtrsim
10^{14}\,M_{\odot}$ at $z=0$ (Toshikawa:2018). Second, any real detection of
protoclusters will depend on their properties at the redshift they are
observed, not on their properties at $z=0$. For instance, a cluster with mass
below our threshold at $z=0$ might have evolved from a very massive
protocluster that happened to have a sustained period of no accretion. While
we include some of the effects of different evolutionary histories (i.e. by
averaging across an ensemble of progenitor objects), our treatment is not
fully realistic because it ultimately relies on a cluster selection at $z=0$.
Finally, obtaining CMB observations of protoclusters detected in a galaxy
survey will require overlapping galaxy and CMB survey footprints. Fortunately,
future galaxy surveys including the Vera C. Rubin Observatory Legacy Survey of
Space and Time (LSST; LSST), Euclid (Euclid) — will have significant overlap
with CMB-S4, and have the potential to detect protoclusters (Brinch:2023;
Araya-Araya:2021). Since the focus of this work is on using CMB surveys to
detect and characterize protoclusters, we postpone a more detailed treatment
of protocluster detection with galaxy surveys to future work. Given some
estimated number of protoclusters detected by a galaxy survey, one can simply
scale our estimated errorbars by $\sqrt{324/N_{c}}$, where $N_{c}$ is the
desired number of protoclusters.

##  3 Results

Figure 7: The fraction of protoclusters expected to have an SZ signal with
${\rm SNR}>5$ at the indicated redshifts as a function of halo mass at $z=0$.

(a) Sunyaev Zel’dovich effect (b) CMB lensing

Figure 8: Stacked SZ (left) and lensing (right) profile measurements. We show
errorbars corresponding to both S4-Wide (dotted) and S4-Ultra deep (solid);
the errorbars for S4-Ultra deep are larger than for S4-Wide because of the
smaller volume probed.

We show the redshift evolution of the SZ and CMB lensing convergence maps in
the absence of noise (but including the instrumental beam) for three different
clusters in Fig. 5. Clusters 1 and 2 (columns one through four) are fairly
typical, while cluster 3 (columns five and six) is the cluster with largest
amplitude SZ signal of all The300 clusters. At high redshift, the SZ and
lensing signals are weak because halos within the protocluster have not
assembled very much mass or gas. By redshift of $z\sim 1$, the three
structures have become approximately virialized galaxy clusters. One can also
see the filamentary structure surrounding the clusters, particularly visible
in the lensing maps at low redshift.

###  3.1 Signal-to-noise of individual cluster measurements

To assess the ability of future CMB surveys to detect protoclusters and to
measure their individual properties, we compute the signal-to-noise of the
mock SZ and lensing profile measurements for each cluster. We define the
signal-to-noise of the $i$th cluster as

${\rm SNR}_{i}\equiv\sqrt{d_{i}^{t}\mathbf{C}\vec{d}_{i}-N_{d}},$ (5)

where $\vec{d}_{i}$ is the profile measurement for the $i$th cluster and
$\mathbf{C}$ is the covariance of the measurement, including both noise and
intrinsic scatter), and $N_{d}$ is the number of elements in $\vec{d}$. We
subtract $N_{d}$ to prevent the estimated SNR from growing arbitrarily large
as the length of the data vector increases, even in the absence of signal.
Note that we do not account for a “look elsewhere” effect associated with
attempting to detect protoclusters at various positions across the sky;
rather, the computed values correspond to the signal-to-noise at the true
cluster location. Note also that the we do not assume a profile shape when
computing the signal-to-noise; the signal-to-noise could in principle be
enhanced with the use of a matched filter.

Fig. 6 shows the distribution of cluster signal-to-noise values for the SZ
(left) and lensing (right) measurements for S4-Ultra deep We have excluded the
highest signal-to-noise cluster (cluster 3 in Fig. 5) when making this figure.
For comparison, we also show the signal-to-noise distribution when the
profiles are computed using only the noise maps (black dashed curves). For
both the SZ and lensing measurements, it is clear that detecting protoclusters
and characterizing their individual properties will be difficult, even with
S4-Ultra deep observations: most protoclusters will not be detected with
signal-to-noise greater than one. For the SZ measurements, which generally
have higher signal-to-noise than the lensing measurements, there is a clear
trend in the SNR distribution towards higher values as one goes to lower
redshift, reflecting the build up of gas pressure in the main halos of the
protoclusters over time. [Note that the evolution with redshift of the SZ
signal-to-noise ratios is faster than for the lensing signal, since $Y\sim
M^{5/3}$ while $\kappa\sim M$. The steep dependence of the SZ signal with halo
mass results from its dependence on both the electron number density, $n_{e}$,
and temperature, $T$, which both increase with increasing halo mass ($Y\sim
n_{e}T$, $n_{e}\sim M$ and $T\sim M^{2/3}$, see Kaiser:1986).]

Although most protoclusters will not be detectable at high significance with
S4-Ultra deep, some structures with large SZ signals may be detected: In Fig.
6, a small but non-zero fraction of protoclusters will be detected with
signal-to-noise greater than five. To make this clearer, Fig. 7 shows the
fraction of protoclusters that will have high-significance SZ detections
(${\rm SNR}>5$) as a function of their mass at $z=0$. The significance of the
SZ detection in this figure depends on both the noise level in the
measurements (instrumental and foreground), and the evolutionary histories of
the protoclusters. Protoclusters that assemble a large fraction of mass at
early times, for instance, will be easier to detect at high redshift. We find
that for clusters with $M\gtrsim 10^{15}\,M_{\odot}/h$ at $z=0$, the chances
of obtaining a high-significance detection of the progenitor structure at
$z\sim 2$ are roughly 20%. We find that essentially no protoclusters will have
high-significance CMB lensing detections. This finding appears to be roughly
consistent with Madhavacheril:2017 and raghunathan22b. The typical mass of the
main progenitor halos in our protocluster sample at $z\sim 2$ is
$10^{13.6}\,M_{\odot}$; both Madhavacheril:2017 and raghunathan22b predict
that clusters with such masses will be below the detection threshold of
CMB-S4. Note, though, that both of those works assume models for the cluster
SZ signal that are based on virialized low-redshift clusters.

###  3.2 Stacked profile forecasts

Since we have shown above that it will be difficult to detect and characterize
individual protoclusters with future CMB observations, we now consider the
prospects for a stacking analysis in which the protocluster measurements are
averaged across many structures to improve signal-to-noise.

In Fig. 8, we show the stacked lensing (right) and SZ (left) profiles
measurements for the protoclusters along with forecast errors for S4-Wide and
S4-Ultra deep. For S4-Wide, we forecast stacked SNR for the SZ measurements of
$[16.5,13,7.2,3.2]$ at redshifts of $[2.0,2.5,3.0,3.5]$, respectively. [We
find that the stacked measurements from S4-Wide will have somewhat higher
precision than those of S4-Ultra deep given the larger number of clusters
probed by S4-Wide (see discussion in §2.5), although we emphasize that this
result is sensitive the assumptions made in §2.5.] For the stacked lensing
measurements, the corresponding SNR values are $[11,7.3,5.7,4.6]$. The signal-
to-noise values increase with decreasing redshift as the main halos of the
protoclusters accrete more mass and gas. We emphasize that it is these main
halos to which the lensing and SZ signals are primarily sensitive: signals
from mass outside these main halos makes a small contribution to the SZ and
lensing signals, and in any case, this contribution is largely removed via the
aperture subtraction process. [As explained previously, our forecasts assume
statistically isotropic noise across the sky, and for this reason, we weight
all clusters equally when computing the stacked profiles. In principle,
signals associated with the protoclusters themselves could lead to variation
in the effective noise at the location of each protocluster. However, since
the individual protocluster signals are essentially undetectable, this is
expected to be a small effect.]

Interestingly, the lensing SNR in the highest redshift bin is actually higher
than that for the SZ because of the rapid evolution of the SZ signal with halo
mass noted previously. The stacked lensing signal can also be detected to
larger radii than the SZ signal for essentially the same reason: since the
number density and temperature of the electron gas both decline with cluster-
centric radius, the SZ signal declines faster with radius than the lensing
signal. The fact that the CMB lensing signal can be detected out to roughly
2.5 Mpc reflects the extended nature of protoclusters.

Figure 9: Relation between the true mass of the main halos within
protoclusters, $M$, and the projected mass $M^{2D}$ that determines the
gravitational lensing signal. At high redshift, when protoclusters are still
very extended, the scatter between 3D mass and projected mass becomes larger.

Figure 10: Comparison of measured $Y$-$M$ relation to predictions based on the
model of Baxter:2023 (B23), which includes the prescription for nonthermal
pressure support from Shaw:2010 (S10) (orange curves in both panels). Left
panel shows the evolution of $Y$ with redshift for individual (main) halos
from The300 (grey curves) and the average across all (main) halos in the
The300 (blue points). The right panel shows the same, but normalized to a
model (black dashed curve) that fixes the nonthermal pressure fraction to its
value at $z=0$. The green band represents the forecast uncertainties on $Y$
from S4-Wide.

Future SZ and CMB lensing protocluster profile measurements can be used to
constrain the integrated $y$ signal, $Y$ (Eq. 2) and main halo mass, $M$,
respectively. We emphasize again that we are primarily concerned here with the
properties of the main halo in the protocluster, which itself may be composed
of many halos. We compute $Y$ within an aperture of radius $R_{\rm
max}=R_{200c}$. While $Y$ can be directly measured from the SZ observations,
the lensing convergence measurements are sensitive to the line-of-sight
projection of the protocluster mass, $M^{2D}$, rather than the true 3D mass:

$M^{2D}=2\pi\int_{0}^{R_{\rm max}}\int dz\,r\,dr\,\rho(r,z),$ (6)

where $\rho(r,z)$ is the density at position $z$ and along the ray to the
cluster, and $r$ is the radial separation from that ray. As with $Y$, We
choose an aperture size of $R_{\rm max}=R_{200c}$. For low-redshift,
virialized clusters, we expect $M^{2D}$ and $M$ to be tightly correlated,
given that clusters are roughly spherically symmetric. Scatter in the relation
between $M$ and $M^{2D}$ is expected due to e.g. triaxiality (Rasia:2012).
High-redshift protoclusters, on the other hand, may be far from spherical. In
Fig. 9 we show a comparison between $M$ and $M^{2D}$ as a function of redshift
[for all simulated clusters]. The 2D mass is always larger than the 3D mass,
since the volume probed is larger by definition. As one goes to higher
redshift, the 2D masses exhibit more scatter around $M$, reflecting the large
asphericity of the protoclusters at high redshift. The scatter around the
1-to-1 relation is roughly $\sigma(\log_{10}M^{2D})=[0.05,0.11,0.12]$ at
$z=[0,2,5]$, respectively. In all cases, this scatter is lower than that due
to noise in the CMB lensing measurements. In generating this figure, we
exclude several clusters for which their evolution leads to them being near
the edge of the high-resolution simulated region at high redshift.

###  3.3 Constraining nonthermal pressure support and virialization

Mergers and accretion events can drive significant bulk motion and turbulence
in the cluster gas. Bulk gas motions can in turn provide significant,
nonthermal pressure support to the gas. Unlike thermal motions, these bulk
motions do not contribute to the thermal SZ signal. Consequently,
understanding the level of nonthermal pressure support in clusters, as well as
its redshift evolution, is critical for many cosmological analyses of galaxy
clusters (e.g Eckert:2019). Measuring the SZ signals of protoclusters and
comparing to theoretical models therefore provides an important test of our
understanding of cluster virialization, nonthermal pressure support, and the
extent to which clusters depart from self-similar expectations at high
redshift.

Fig. 10 shows the evolution of $Y$ as a function of redshift for the main
progenitor halos of The300 clusters. The grey curves represent individual
halos, while the blue points represent the average across all halos. It is
clear that the evolution of $Y$ with redshift varies significantly between
different clusters owing to their different evolutionary histories. The orange
curve is the model which we describe below.

To model the evolution of $Y$ with redshift, we adopt the model for the
$Y$-$M$ relation from Baxter:2023, which in turn draws heavily from
KomatsuSeljak2001; Komatsu:2002; we refer readers to those works for more
details of the modeling choices, and provide only a summary here. The
KomatsuSeljak2001; Komatsu:2002 model assumes that mass distribution of the
halo is described by a Navarro-Frenk-White profile (NFW), and that the halo
gas is in hydrostatic equilibrium within this potential. The amount of gas in
the halo is set equal to the cosmic baryon fraction, and the gas density
profile is assumed to follow the matter profile at large cluster-centric
distances. The Baxter:2023 model additionally includes a prescription for
decline in the halo pressure profile due to virial shocks at the halo
splashback radius (Bertschinger:1985; Diemer:2014), which can have a
significant impact on the $Y$-$M$ relation, and also a simple prescription for
the two-halo term, which accounts for the contributions to $Y$ from nearby
halos (e.g. Vikram:2017). Most importantly for our purposes, the Baxter:2023
model includes a prescription for nonthermal pressure support which is in turn
derived from Shaw:2010. Briefly, the fraction of the total pressure supplied
by nonthermal sources is assumed to have a radial and redshift dependence of
the form

$f_{\rm nth}(r,z)=\alpha_{\rm nt}(z)(r/R_{500c})^{0.8},$ (7)

where

$\displaystyle\alpha_{\rm nt}(z)$ $\displaystyle=$
$\displaystyle\alpha_{0}f(z)$ (8) $\displaystyle f(z)$ $\displaystyle=$
$\displaystyle{\rm min}\left[(1+z)^{\beta},(f_{\rm max}-1){\rm tanh}(\beta
z)+1\right]$ (9) $\displaystyle f_{\rm max}$ $\displaystyle=$ $\displaystyle
4^{-n_{\rm nt}/\alpha_{0}},$ (10)

with $\alpha_{0}=0.18$ and $n_{nt}=0.8$. This model was found by Shaw:2010 to
provide a good match to the simulations of Nagai:2007. Shaw:2010 also found
weak evidence for redshift evolution of $n_{\rm nt}$ as well as dependence of
the redshift evolution of nonthermal pressure on the details of the baryonic
modeling.

The orange curves in both panels of Fig. 10 represent the calculation of $Y$
using the Baxter:2023 model with the Shaw:2010 prescription for nonthermal
pressure support. In calculating this model curve, we adopt the true mean mass
of the main halos of the protocluster structures at each redshift. We find
that this model provides a remarkably accurate description of the average
redshift evolution of $Y$ over a wide range of redshift (left panel of Fig.
10). This is perhaps surprising given that protoclusters can significantly
violate many of the assumptions of the Baxter:2023 model. It appears, though,
that while individual protoclusters may exhibit significant departures from
e.g. a NFW density profile leading to significant scatter in their $Y$
evolution, these differences tend to wash out when averaged over many
clusters. In the right panel of Fig. 10 we plot the redshift evolution of our
baseline model (again with the orange curve) relative to a model for which we
fix $\alpha_{\rm nt}(z)=\alpha_{0}f(0)$, i.e. for which there is no redshift
evolution in the nonthermal pressure fraction. Between redshifts of $z\sim 2$
to $z\sim 5$, redshift evolution of the nonthermal pressure fraction is
expected to suppress the mean $Y$ by about 10%. This is because in the model
of Shaw:2010 (see also Stanek:2010), the nonthermal pressure fraction
increases with redshift (up to some maximum value), leading to a suppression
of $Y$. The green band in the right panel of Fig. 10 shows the expected
uncertainty on $Y$ given the stacked measurements from S4-Wide. It appears
that the precision offered by CMB-S4 will be sufficient to test a plausible
range of predictions for the redshift evolution of nonthermal pressure
support. Note, though, that uncertainties on mass estimates for protoclusters
may complicate such an analysis. We postpone a detailed investigation of this
possibility to future work. The blue curve in the right panel of Fig. 10
represents the evolution of $Y$ as measured in The300. There is significant
scatter in these measurements, particularly at high redshift, due to intrinsic
variation between the clusters (see also the left panel of the figure). To
within this scatter, The300 clusters agree reasonably well with the redshift
evolution of $Y$ predicted by the Shaw:2010 model, although there may be a
hint of slightly less redshift evolution of the nonthermal pressure support.

##  4 Summary

We have explored the ability of future CMB surveys like CMB-S4 to detect and
characterize the properties of high-redshift protoclusters. Relative to
previous studies, a unique feature of our analysis is that we use realistic
simulated protoclusters from The300 in our forecasts, rather than analytic
models based on low-redshift clusters. This is important because the
assumptions of e.g. spherical symmetry and hydrostatic equilibrium can be
grossly violated for protoclusters.

We find that it will be challenging to detect the lensing and SZ signals from
individual protoclusters, even with high-sensitivity data from S4-Ultra deep.
However, we do expect that of order 10% of the progenitors of very massive
clusters ($M\gtrsim 10^{15}M_{\odot}$) will have SZ signals detectable at high
significance at $z=2.5$ (see Fig. 7). It is unlikely that any protoclusters at
high redshift will be detected via their CMB lensing signals.

To date, most protoclusters have been detected via overdensities of high-
redshift galaxies. Assuming that surveys like Euclid, LSST and the Nancy Grace
Roman Space telescope (Spergel:2015) can identify samples of protoclusters via
galaxy overdensities in the region surveyed by CMB-S4, it will be possible to
constrain their average gas content and mass using CMB-S4 (Fig. 8). Assuming
all protoclusters that are progenitors of clusters at $z=0$ with $M\gtrsim
10^{15}\,M_{\odot}$ can be identified, we find that the CMB lensing signal,
which is directly related to the cluster mass, can be measured with a signal-
to-noise of roughly 5.7 at $z\sim 3$, while the SZ signal can be measured with
signal-to-noise of 7.2 at this redshift. Because of their non-spherical
nature, protoclusters will exhibit significantly more scatter in their lensing
(i.e. projected) masses than low-redshift virialized clusters (Fig. 9).
However, this scatter will be small compared to the uncertainties from the CMB
lensing mass estimates.

Finally, the measurements of the protocluster SZ signals with S4-Wide can be
used to test models for the redshift evolution of nonthermal pressure support
within protoclusters (Fig. 10), providing insight into the processes by which
galaxy clusters assemble and virialize.

We emphasize that several of the forecasts in this work are likely
conservative, since our stacking analysis is predicated on being able to
identify only those protoclusters which are the progenitors of clusters at
$z=0$ with $M\gtrsim 10^{15}\,M_{\odot}$. It is likely that future wide-field
galaxy surveys with significant overlap with S4-Wide will be able to detect
less massive structures. Our forecasts can be trivially scaled to account for
any desired number of protoclusters; we postpone investigation of more
realistic protocluster selections to future work.

##  5 Data availability

This manuscript was developed using data from The300 galaxy clusters sample.
The data is available on request following the guideline of The300
collaboration (https://www.the300-project.org). The data used to make the
figures shown in this work are available upon request.

## Acknowledgments

This work was partially supported by the Center for AstroPhysical Surveys
(CAPS) at the National Center for Supercomputing Applications (NCSA),
University of Illinois Urbana-Champaign.

This work made use of the Illinois Campus Cluster, a computing resource that
is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with
the National Center for Supercomputing Applications (NCSA) and which is
supported by funds from the University of Illinois at Urbana-Champaign. This
work also used the computational and storage services associated with the
Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and
Education’s Research Technology Group.

We thank Lukas Zalesky for fruitful discussions related to this work.

WC is supported by the STFC AGP Grant ST/V000594/1 and the Atracción de
Talento Contract no. 2020-T1/TIC-19882 granted by the Comunidad de Madrid in
Spain. He also thanks the Ministerio de Ciencia e Innovación (Spain) for
financial support under Project grant PID2021-122603NB-C21 and ERC: HORIZON-
TMA-MSCA-SE for supporting the LACEGAL-III project with grant number
101086388.


# Branched flows of flexural elastic waves in non-uniform cylindrical shells

Kevin Jose  Faculty of Engineering and Physical Sciences, University of
Southampton, Southampton, SO17 1BJ, UK  Neil Ferguson  Faculty of Engineering
and Physical Sciences, University of Southampton, Southampton, SO17 1BJ, UK
Atul Bhaskar  Faculty of Engineering and Physical Sciences, University of
Southampton, Southampton, SO17 1BJ, UK  Department of Mechanical Engineering,
University of Sheffield, Mappin Street, Sheffield, S1 3JD, UK
<EMAIL_ADDRESS>
###### Abstract

Propagation of elastic waves along the axis of cylindrical shells is of great
current interest due to their ubiquitous presence and technological
importance. Geometric imperfections and spatial variations of properties are
inevitable in such structures. Here we report the existence of branched flows
of flexural waves in such waveguides. The location of high amplitude motion,
away from the launch location, scales as a power law with respect to the
variance and linearly with respect to the correlation length of the spatial
variation in the bending stiffness. These scaling laws are then theoretically
derived from the ray equations. Numerical integration of the ray equations
also exhibit this behaviour—consistent with finite element numerical
simulations as well as the theoretically derived scaling. There appears to be
a universality for the exponents in the scaling with respect to similar
observations in the past for other types of waves, as well as flexural and
hence dispersive waves in elastic plates.

## Introduction

Waves propagating through heterogeneous media with spatially correlated
randomness show a peculiar behaviour, known as branched flows[1], in a variety
of physical contexts such as optics[2], microwaves[3], electron waves[4],
tsunami waves[5], sound waves in ocean[6] etc. The phenomenon of branched
flows of waves is characterised by the emergence of flow-like patterns with
spatial branching. Emergence of branching is also associated with the
occurrence of focusing events which leads to regions of high amplitude.
Additionally, the expected distance $\langle l_{f}\rangle$ – of the first of
such focusing events from the point of launch – scales as

$\langle l_{f}\rangle\propto L_{\mathrm{c}}\langle h^{2}\rangle^{-1/3},$ (1)

where $\langle\cdot\rangle$ signifies the mean, $L_{\mathrm{c}}$ is the
correlation length of the isotropic randomness field and $\langle
h^{2}\rangle$ is a non-dimensional measure of the severity of the randomness;
$h$ is defined more rigorously below.

In the context of elastic waves, we recently showed the existence of branched
flows and the associated scaling law in thin elastic plates [7]. This raises a
natural question about the universality of branched flows in mechanical waves
carried by elastic structures. In many applications such as health monitoring
of buried gas pipes, shells act as waveguides for elastic waves, hence they
provide a practical context. Cylindrical shells also provide a practical
context of an elastic waveguide to study branched flows without reflections
from the edges [7]. Further, the need for contrived periodic boundary
conditions, as previously implemented [8], is obviated. Waves in elastic
shells have attracted the attention of dynamicists for sometime [9, 10].
Pioneering work on the statics of cylindrical shells was carried out by
Donnell[11]. This has been extended, for example, by Yu[12, 13] and Naghdi &
Cooper[14], to derive the complete equations of motion. Numerous formulations
of the dynamics of shells, of varying degrees of accuracy, exist; these have
been summarized by Greenspon[10] and Leissa[15]. Propagating waves[16, 17,
18], normal modes[12, 13], shells under random excitation [19] have been
studied theoretically, computationally[20] and experimentally[21]. The
propagation behaviour of plane waves in the presence of inevitable
manufacturing tolerances has not been studied in any detail so far. Here we
examine the the effect of such spatial non-homogeneities and explore the
emergence of channels of energy flow in such elastic waveguides.

We consider the propagation of flexural waves through a hollow cylindrical
waveguide when the wavelength of interest is much shorter than the correlation
length of the heterogeneity of properties, e.g. thickness or material
stiffness. The tubular cross-section undergoes breathing displacement that is
radial, as a deformation pattern propagates axially as a wave. The nominal
thickness of the elastic cylinder is much smaller than the wavelength of
interest, in order to justify ignoring shear through the thickness in our
analyses. Consider a cylinder of non-uniform thickness with the axis along the
$x$-direction; the circumferential direction along $s$ (with units of length)
and thickness $H(x,s)$ about the nominal cylindrical surface. The thickness of
the hollow elastic cylinder has the form $H(x,s)=H_{0}(1-h(x,s))$ where
$h(x,s)$ is a smoothly varying random field with an isotropic auto-correlation
function; the correlation length is $L_{\mathrm{c}}$, and $\langle
h\rangle=0,\ \langle h^{2}\rangle\ll 1$ (see Figure 1).

Figure 1: Spatial variation of thickness over the surface of a cylindrical
waveguide. The radius of the cylindrical shell is $R$, axial length $L_{x}$,
thickness $H$ that varies spatially as per $H(x,s)=H_{0}(1-h(x,s))$, where
$\langle h\rangle=0,\langle h^{2}\rangle\ll 1$. Heatmap shows the spatial
variation of $h(x,s)$. Inset on the top left is an “unwrapped” view of
$h(x,s)$. In the middle of the inset, a circle of radius $L_{\mathrm{c}}$,
equal to the correlation length, is shown.

## Results

Consider an initially plane wavefront of predominantly one wavelength
propagating axially through a thin cylindrical shell. The assumption of slow
spatial variation of properties enables the simplification of the wave
elastodynamics to a set of ray equations using the eikonal/WKB
approximations[22, 23]. They are further simplified using the paraxial
approximation[24], permitted by the weak scattering nature of the problem,
which asserts a predominantly axial direction of the wave vector. These ray
equations are then used to analytically derive the scaling law (Equation 1)
relating the position of focusing and the severity of the non-homogeneity. We
also numerically integrate the ray equations to investigate the emergence of
branched flow and to validate the theoretically derived scaling law. The same
scaling is also probed using finite element simulations which capture the
complete wave elastodynamics, including dispersion.

### Ray equations

#### Thin shell theory

Plane waves of predominantly single wavelength $k_{0}$ are launched at the
left end of the elastic waveguide. Rays are fully described by four
quantities: spatial variables $x,s$ indicating the location of the ray along
the axis and circumference respectively, and $k_{x},k_{s}$, the wavenumber
components in the axial and circumferential directions respectively. The
spatial variables are non-dimensionalised with respect to the nominal radius,
i.e. $\tilde{x}=x/R,\tilde{s}=s/R$, whereas the wavenumbers are non-
dimensionalised with respect to the initial wavenumber
$\tilde{k}_{x}=k_{x}/k_{0},\ \tilde{k}_{s}=k_{s}/k_{0}$. Given the parameter
regime of interest, (thin shell, small curvature, and slowly varying
parameters) we can use the dispersion relation given by Pierce [16] to obtain
the non-dimensionalized ray equations

$\partial_{\tau}\tilde{x}=\tilde{k}_{x}\left(\tilde{k}_{x}^{2}+\tilde{k}_{s}^{2}+\alpha\frac{(-\tilde{k}_{x}^{2}\tilde{k}_{s}^{2})}{\left(\tilde{k_{s}}^{2}+\tilde{k_{x}}^{2}\right)^{3}}\right),\quad\partial_{\tau}\tilde{s}=\tilde{k}_{s}\left(\tilde{k}_{x}^{2}+\tilde{k}_{s}^{2}+\alpha\frac{(\tilde{k}_{x}^{4})}{\left(\tilde{k_{s}}^{2}+\tilde{k_{x}}^{2}\right)^{3}}\right),\quad\begin{bmatrix}\partial_{\tau}\tilde{k}_{x}\\\
\partial_{\tau}\tilde{k}_{s}\end{bmatrix}=\frac{1}{2}(\tilde{k}_{x}^{2}+\tilde{k}_{s}^{2})^{2}\tilde{\nabla}\tilde{h},$
(2)

where $\alpha=12(\nu^{2}-1)\frac{\gamma^{2}}{\epsilon^{2}},\
\gamma=(k_{0}R)^{-1},\ \epsilon=H_{0}k_{0}$. Additionally,
$\tilde{h}(x/R,s/R)=h(x,s)$ and $\tilde{\nabla}=[\partial_{\tilde{x}}\
\partial_{\tilde{s}}]^{T}$. $\tau$ is an arbitrary time scale, consistent with
ray approximations. The initial condition is $\tilde{x}=0,\ \tilde{k}_{x}=1,\
\tilde{k}_{s}=0$ and $\tilde{k}_{s}\in[0,2\pi R)$ depending on circumferential
location of the ray’s starting point.

#### Paraxial ray equations

In the regime of weak scattering studied here, the thin shell ray equations
can be simplified. $\tilde{k}_{x}$ and $\tilde{k}_{s}$, the wavenumbers in the
axial and circumferential directions respectively, are not expected to vary
substantially from their initial values of 1 and 0 respectively, as the
initially launched plane wave is purely axial. For weak scattering of
initially plane waves, it is customary to make the simplifying assumption that
the wave-number does not change at all in the main propagation direction. This
can be achieved by setting $\partial_{\tau}\tilde{k}_{x}=0,\ \tilde{k}_{x}=1$
in Equation 2, which is the essence of the paraxial approximation. Finally,
dropping all $\tilde{k}_{s}$ terms compared to $\mathcal{O}(1)$ terms, the ray
equations become

$\partial_{\tau}\tilde{x}=1,\quad\partial_{\tau}\tilde{s}=\tilde{k_{s}}(1+\alpha),\quad\partial_{\tau}\tilde{k}_{x}=0,\quad\partial_{\tau}\tilde{k}_{s}=\frac{1}{2}\partial_{\tilde{s}}\tilde{h}.$
(3)

The approximations made to obtain Equation 3 from the full thin shell ray
equation use arguments about the physics of the problem. Nonetheless, the
validity of these assumptions is further confirmed later here by numerical
integration of ray equations (details in Methods).

Figure 2: Comparison of thin shell ray equations and equations obtained from
making the paraxial approximation. Left: Ray propagation using full thin shell
theory and paraxial approximation of thin shell theory. Circular and cross
markers indicate location of the first caustic. The $\tilde{s}$ axis has been
“unwrapped” for representational purposes. Middle: Same ray propagation shown
on the cylindrical shell. Right: Plot of the temporal evolution of quantities
describing the one of the rays. Here, markers indicate temporal location of
the first caustic. It can be seen that, at higher values of $\langle
h^{2}\rangle$, the location of the first caustic detected from the paraxial
approximation differs from thin shell.

Figure 2 (left) shows a comparison between the rays obtained from numerical
integration of the full thin shell ray equations (Equation 2), and the
simplified equations obtained following the paraxial approximation (Equation
3) for increasing values of $\langle h^{2}\rangle$. The rays have been plotted
“unwrapped” in the left column of the figure, where 300 rays equi-spaced along
the circumference at the left end are launched as they curve and veer
downstream due to scattering. Transmission behavior on the cylindrical
surface, is shown in Figure 2 (middle) where interesting spiral structures
with no preferred handedness, as expected, are observed. The rays as computed
using thin shells theory versus that using the paraxial approximation look
fairly similar, confirming the validity of the paraxial approximation. It can
be seen that, at the time instant at which the simulation is terminated, all
rays have reached the right edge in case of the paraxial approximation, unlike
the rays obtained from integrating the thin shell ray equations. This is a
consequence of assuming $\partial_{\tau}\tilde{k}_{x}=0$ under the paraxial
assumption. Most importantly, it can be seen that the spatial location of the
first caustic, indicated by circular and cross-shaped markers, is
approximately the same from both sets of ray equations, except for the highest
levels of $\langle h^{2}\rangle$ shown here.

In Figure 2 (right), the time evolution of the four quantities describing one
ray is plotted. We can see that the values of $\tilde{x},\ \tilde{s},\
\tilde{k}_{s}$ are in excellent agreement. $\tilde{k}_{x}$ has a constant
value in the paraxial approximation as stated earlier. While it does not show
the variation with time that the thin shell case shows, note that the values
do not change appreciably from $1$. Regardless, the first caustic is detected
by looking for the instant where the $\tilde{s}-\tilde{k}_{s}$ curve becomes
locally two valued (see Methods) and is not dependant on $\tilde{k}_{x}$
directly. In Figure 2 (right), markers indicate temporal location of the first
caustic.

The paraxial approximation becomes more inaccurate as time passes. However,
since we are interested only in the location of the first caustic, and they
tend to appear fairly early, this eventual drift is inconsequential. The
paraxial approximation also breaks down faster when the value of $\langle
h^{2}\rangle$ is higher, i.e. when the weak scattering assumption starts
breaking down. However, this is partially counteracted by the fact that when
the value of $\langle h^{2}\rangle$ is higher, the first caustic appear
earlier. Nonetheless, the progressive degradation of the paraxial
approximation with increase in $\langle h^{2}\rangle$ can be seen in Figure 2.

We will now use these ray equations to derive the scaling of $\langle
l_{f}\rangle$ analytically. We will also numerically integrate them to study
the scaling. The ray equations Equation 2 (and, therefore, Equation 3 too) are
obtained from simplified dispersion relation for flexural waves in a thin
shell. At the higher $\langle h^{2}\rangle$, some rays are expected to
completely back-scatter due to the severity of randomness. This is not
captured by Equation 2. In the Supplementary Information, we derive the ray
equations starting from the governing equations of the displacements of a thin
cylindrical shell. These ray equations are more sophisticated and show the
expected back scattering at higher values of $\langle h^{2}\rangle$. However,
after applying the paraxial approximation, the resultant set of ray equations
show very similar results to Equation 3. The analytical and numerical
validation of the scaling law Equation 1 using these ray equations is also
shown in the Supplementary Information.

### Scaling law: from analysis, ray equations, and finite element
elastodynamics

We use the ray equations obtained after making the paraxial approximations
Equation 3, and following a process very similar to other branched flow
works[7, 5] we obtain, $\langle l_{f}\rangle\sim L_{c}\langle
h^{2}\rangle^{-1/3}$ barring an extremely weak dependence of the
proportionality constant on $\gamma$ and $\epsilon$ (see Supplementary
Information).

The scaling of the location of the first caustic resulting from the analysis
of ray equations can also be obtained numerically from: (i) numerical
integration of the ray equations, and (ii) finite element elastodynamic
simulations. The emergence of branched flows is clearly visible from finite
element elastodynamics simulations; an example of which is shown in Figure 3.
‘Snapshots’ of the temporal evolution are shown in (a $\to$ b $\to$ c $\to$ d)
of an initially plane wave front, as it propagates along the elastic cylinder
with non-uniform thickness. The entire domain is shown on the top of each
panel and regions of high amplitude, indicated by colored lines, are zoomed
into and shown on the bottom of each panel. The radial displacement has been
greatly exaggerated for representational purposes. Note that the absolute
value of displacement can be scaled arbitrarily since we are considering the
linear regime. The initially plane wave front (a) splits into distinct
branches (b) leading to regions of extreme amplitudes. As the wavefront
propagates further (c, d) more branching is observed. The widening of the
wavefront as expected because of the dispersive nature of flexural waves in
shells as seen here.

Figure 3: Emergence of branched flow in for an initially plane wavefront in a
thin cylinder of non-uniform thickness. Temporal evolution (a $\to$ b $\to$ c
$\to$ d) of an initially plane flexural wave front in a thin elastic cylinder
of non-uniform thickness. The full cylinder has been shown on the top of each
panel. The regions of high amplitude at each time instant is indicated with
colored lines and they been zoomed into and shown at the bottom of each panel.
The emergence of branching leading to locations of extreme amplitudes and
widening of the wavefront consistent with the dispersive character of this
class of elastic waves is also observed.

Using the approach of detecting the first caustic described in the Methods
section, the scaling of location of the first caustic with the statistical
properties of the random field and the geometry of the cylinder are now
explored. Figure 4 shows results for the distance of the first focusing event
from the point of launch as a function of the severity of randomness for 800
realisations of the shells of correlated randomness. The mean for each value
of $\langle h^{2}\rangle$ is shown by a blue square marker. The figure of the
left is from the numerical integration of ray equations whereas that on the
right is from FE simulations, both showing good agreement with the scaling
$\langle l_{f}\rangle\sim\langle h^{2}\rangle^{-1/3}$.

Figure 4: Scaling of location of the expected location of the first focusing
event with “severity” of randomness. Locations of first caustic as obtained
from numerical ray integration (left) and FE elastodynamics simulations
(right). Circular markers indicate locations of first caustic from individual
simulations. Square markers show the expected location of the first caustic.
Vertical lines indicate 2 standard deviation (centered around the mean). The
expected locations of the first caustic show the expected power law scaling
with $\langle h^{2}\rangle$.

The linear scaling of the location of focusing with the correlation length
$\langle l_{f}\rangle\propto L_{\mathrm{c}}$ is also validated using the two
numerical methods, see Figure 5. Using the curve obtained for
$L_{\mathrm{c}}/L_{\mathrm{c}}^{\mathrm{ref}}=1$ as reference, the prediction
assuming the linear scaling with $L_{\mathrm{c}}$ is shown using dashed lines.
It is clear that the actual simulations (markers) agree with this prediction
both for ray simulations as well as finite elements elastodynamics.. In the
results from FE elastodynamics simulations, the scaling seems to diverge from
the prediction at higher values of $\langle h^{2}\rangle$; this is expected
since the weak scattering assumption breaks down at these values of $\langle
h^{2}\rangle$. This does not happen in the simulations using numerical
integration of ray equations since we use the formulation obtained from
applying the paraxial approximation which assumes weak scattering. The linear
scaling can also be inferred from dimensional arguments that the only length
scale in the problem is the correlation length.

Figure 5: Scaling of the expected location of the first focusing event with
correlation length of the randomness. Expected locations of the first caustic
from numerical ray integration (left) and FE elastodynamics simulations
(right) for different correlation lengths. Taking the curve corresponding to
$L_{c}/L_{c}^{\mathrm{ref}}=1$ to be the reference, dashed lines show the
predicted behavior at other correlation lengths assuming $\langle
l_{f}\rangle\propto L_{c}^{1}$. The actual simulations agree well this
prediction hence confirming the linear scaling.

The insensitivity of $\langle l_{f}\rangle$ to wavelength (for $\lambda\ll
L_{\mathrm{c}}$) is confirmed by the two numerical approaches, see Figure 6.
The expected breakdown of the scaling due to strong scattering is seen here
for the FE elastodynamics simulations. Finally, using ray simulations, we
verified that $\langle l_{f}\rangle$ is insensitive to variations of the
radius of the cylinder in the shallow cylinder regime $\lambda\ll R$, see
Figure 7. The same could not be carried out using FE simulations due to
computational constraints. Note that, since $\tilde{l}_{f}$ is non-
dimensionalised with respect to the radius, a nominal radius $R_{0}$ is used
to scale $\tilde{l}_{f}$ appropriately when studying (in Figure 7) the effect
of changing radius.

Figure 6: Scaling of the expected location of the first focusing event
wavelength of the initial wavefront. Expected locations of the first caustic
from numerical ray integration (left) and FE elastodynamics simulations
(right) for different wavelengths. As long as $\lambda\ll L_{c}$,
$\langle\tilde{l}_{f}\rangle$ is independent of wavelength. The power law
scaling seems to break down at higher $\langle h^{2}\rangle$, especially in FE
simulations. This is consistent with the expectation that the scaling holds
only for weak scattering and higher $\langle h^{2}\rangle$ corresponds to
higher scattering.  Figure 7: Scaling of the expected location of the first
focusing event with the radius of the cylinder. Scaling of the expected
locations of the first caustic with radii from numerical ray integration. The
expected location of the first caustic is independent of the radius in the
parameter ranges of interest.

## Conclusions

There appears to be a universality of branched flows in the propagation of
elastic waves through random media. Following our previously reported results
on branched flow in elastic plates[7], here we demonstrate analogous behaviour
for shells with correlated random properties. The ray equations for flexural
waves in a thin elastic cylinder are derived. A paraxial approximation,
permitted by the parameter regimes of interest, is applied to the ray
equations. This is used to analytically demonstrate that the expected location
of the first caustic in shallow cylinders shows the scaling $\langle
l_{f}\rangle\propto L_{\mathrm{c}}\langle h^{2}\rangle^{-1/3}$. This scaling
was then corroborated using numerical integration of ray equations and full FE
elastodynamic simulations.

An immediate extension of this work on cylindrical elastic shells would be to
explore the dependence of the scaling of the first caustic with the radius for
shells with appreciable curvature (i.e small radius). We are unable to do so
in this work since the requirement of $2\pi R\gtrsim L_{\mathrm{c}}$ limits
how small the radius can be. This can be remedied by using anisotropic
randomness which would enable one to reduce the radius by using a smaller
correlation length in the circumferential direction. There may be an elegant
scaling of $\langle l_{f}\rangle$ with radius in this parameter regime. The
existing literature on branched flows in media with spatially anisotropic
randomness[25] can be leveraged.

## Methods

### Detecting first focusing events

For numerical ray simulations, we use a similar method as the one described in
our earlier work[7] to detect focusing events. The first focusing event is
detected by finding the time and location when the
$\tilde{k}_{\tilde{s}}-\tilde{s}$ curve folds over itself[26, 27]. This is
done numerically by tracking the local slope of the curve and detecting
locations when the slope become higher than $\pi/2$. The temporal evolution of
the $\tilde{k}_{\tilde{s}}-\tilde{s}$ curve is plotted in Figure 8. The first
time this curve folds over itself, signalling a caustic, is shown in red. Note
that since, $\tilde{k}_{\tilde{s}}$ has values centered around $0$, a constant
positive value is added to it for representational purposes. This value is not
used when numerically detecting the caustic.

Figure 8: Plot of the temporal evolution of the
$\tilde{k}_{\tilde{s}}-\tilde{s}$ curve obtained from numerical ray
simulations. Since values of $\tilde{k}_{\tilde{s}}$ are centered around zero,
a constant positive offset is added for the purpose of visualization. The
first caustic is detected by finding the location where the
$\tilde{k}_{\tilde{s}}-\tilde{s}$ curve folds over itself (shown in red).

We use a very similar technique as the one used in existing branched flow
literature [7] to detect the location of the first caustic from FE
simulations. For each simulation, the displacement fields are used to
construct the integrated intensity,
$I(x,s)=\int_{0}^{T}\eta^{2}(x,s,t)\mathrm{d}t$. This is in turn used to
construct the scintillation index, $S(x)=\langle I^{2}\rangle_{s}/\langle
I\rangle^{2}_{s}-1$. The location of the first significant peak of $S(x)$
indicates the location of the first caustic (see Figure 9).

Figure 9: Plot of integrated intensity (top) and scintillation index (bottom)
of an exemplar FE simulation. Note that the integrated intensity has been
normalized along by the mean along the circumferential direction for visually
emphasize the locations of extreme amplitudes. The first prominent peak
(circular marker) of the scintillation index ($SI(x)$) curve indicates the
location of the first caustic ($l_{f}$).

### Remarks on numerical ray integration and FE elastodynamics simulations

According to the paraxial approximation, rays travel at a constant speed
($c_{1}$) in the main propagation direction. This enables us to calculate the
total time the rays would take to travel the entire length of the cylinder.
The time stepping in the numerical ray simulations is set using this to ensure
that there are 2000 time steps in the full length traversal. We can be certain
that this time stepping is adequate since no variation is seen in the rays
when the time steps are varied around 2000. Another way of being confident of
the adequacy of the time stepping is to note that the time steps roughly
corresponds the same number of steps the main propagation direction and 2000
steps is adequate to capture the features in the main propagation direction
which are of the order of the correlation length i.e ($L_{x}/2000\ll
L_{\mathrm{c}}$). Note that we are unable to use any other physical arguments
to arrive at the time step since the time scaling for ray equations is
arbitrary as shown earlier.

For the FE elastodynamics simulations, a combination of physical arguments and
some trial and error is used to find the total time of the simulations. The
time step is set to $\frac{\pi}{4\omega}$, this ensures that there are 8
points per cycle since the forcing has predominant frequency component of
$\omega$. The cylinder is meshed using rectangular element which have the
dimensions $\lambda/8$ and $3\lambda/8$ in the axial and circumferential
directions respectively. Note that the meshing can be a bit coarser in the
circumferential direction since the wave predominantly travels in the axial
direction and therefore the variations along the circumference are modest and
of length scales much longer than the wavelength. The element aspect ratio of
1:3 is typically inadvisable in general, however, since we are confident that
the variations in the circumferential direction are modest and of longer
length scales, this aspect ratio will not lead to ill conditioning or mis-
estimation of results. The spatial discretization of $\lambda/8$ ensures that
the spatial variation due to a wave of predominant wavelength $\lambda$ is
adequately captured. Admittedly, since the system being studied is dispersive,
wavelength shorter than $\lambda$ will also be excited and the chosen
discretization may not model them well. The computational expense of these FE
simulations must be emphasized here and this is the reason behind some
compromises in the choices made during meshing.

It was ensured that the location of the detected caustics from numerical ray
integration did not show any angular bias. Figure 10 show the angular
distribution of caustics detected using numerical ray integration. No
consistent and appreciable bias is visible. Angular bias can be introduced
inadvertently from improper interpolation of the randomness field when
conducting numerical ray integration.

Figure 10: Histograms of angular locations of the first focusing events
obtained from numerical ray integration. Polar histograms of the angular
locations, shown in blue, show no consistent appreciable angular bias. The
green circle indicates the outline of a polar histogram with perfectly zero
angular bias. The individual data-points used to construct the histograms are
shown with orange dots. Their angular location is the one obtained from
simulations; their radial position in the above plots is randomized in the
interest of visualization.

Ray integration and FE simulation routines developed for our work on branched
flows in elastic plates[7] were suitably modified for the elastic cylindrical
domain under consideration. Most of the required modifications pertain to
imposing the continuity condition as one went around the circumferential
direction. The method detailed in [7] for generating random field of specified
correlation length, fortuitously, generates fields which are periodic along
the parallel sides and hence, that code was ready to repurposed with minimal
modification for generating $h(x,s)$ since the continuity requirement was
automatically satisfied. The code snippet detecting caustics numerically
during ray simulations was rewritten to respect continuity in the
circumferential direction. The parametrized FE code was modified to generate
the right circular cylindrical geometry. The code to export integrated
intensities was also modified in light of the different coordinate system.

## References

  * [1] Heller, E. J., Fleischmann, R. & Kramer, T.  Branched flow.  _Physics Today_ 74, 44–51 (2021).  https://doi.org/10.1063/PT.3.4902.
  * [2] Patsyk, A., Sivan, U., Segev, M. & Bandres, M. A.  Observation of branched flow of light.  _Nature_ 583, 60–65 (2020).  https://doi.org/10.1038/s41586-020-2376-8.
  * [3] Höhmann, R., Kuhl, U., Stöckmann, H.-J., Kaplan, L. & Heller, E. J.  Freak waves in the linear regime: A microwave study.  _Phys. Rev. Lett._ 104, 093901 (2010).  https://doi.org/10.1103/PhysRevLett.104.093901.
  * [4] Topinka, M. _et al._ Coherent branched flow in a two-dimensional electron gas.  _Nature_ 410, 183–186 (2001).  https://doi.org/10.1038/35065553.
  * [5] Degueldre, H., Metzger, J. J., Geisel, T. & Fleischmann, R.  Random focusing of tsunami waves.  _Nature Physics_ 12, 259–262 (2016).  https://doi.org/10.1038/nphys3557.
  * [6] Wolfson, M. A. & Tomsovic, S.  On the stability of long-range sound propagation through a structured ocean.  _The Journal of the Acoustical Society of America_ 109, 2693–2703 (2001).  https://doi.org/10.1121/1.1362685.
  * [7] Jose, K., Ferguson, N. & Bhaskar, A.  Branched flows of flexural waves in non-uniform elastic plates.  _Communications Physics_ 5, 1–6 (2022).  https://doi.org/10.1038/s42005-022-00917-z.
  * [8] Metzger, J. J., Fleischmann, R. & Geisel, T.  Statistics of extreme waves in random media.  _Physical Review Letters_ 112, 1–5 (2014).  https://doi.org/10.1103/PhysRevLett.112.203903.
  * [9] Soedel, W.  _Vibrations of shells and plates_ (CRC Press, 2004).
  * [10] Greenspon, J. E.  Vibrations of a thick-walled cylindrical shell—comparison of the exact theory with approximate theories.  _The Journal of the Acoustical Society of America_ 32, 571–578 (1960).  https://doi.org/10.1121/1.1908148.
  * [11] Donnell, L. H.  Stability of thin-walled tubes under torsion.  Tech. Rep., National Advisory Committee for Aeronautics/NASA (1935).  https://ntrs.nasa.gov/citations/19930091553.
  * [12] Yu, Y.-Y.  Free Vibrations of Thin Cylindrical Shells Having Finite Lengths With Freely Supported and Clamped Edges.  _Journal of Applied Mechanics_ 22, 547–552 (1955).  https://doi.org/10.1115/1.4011152.
  * [13] Yu, Y.-Y.  Vibrations of thin cylindrical shells analyzed by means of Donnell-type equations.  _Journal of the Aerospace Sciences_ 25, 699–715 (1958).  https://doi.org/10.2514/8.7849.
  * [14] Naghdi, P. M. & Cooper, R. M.  Propagation of elastic waves in cylindrical shells, including the effects of transverse shear and rotatory inertia.  _The Journal of the Acoustical Society of America_ 28, 56–63 (1956).  https://doi.org/10.1121/1.1908222.
  * [15] Leissa, A. W.  Vibration of shells.  Tech. Rep., National Advisory Committee for Aeronautics/NASA (1973).  https://ntrs.nasa.gov/citations/19730018197.
  * [16] Pierce, A. D.  Waves on Fluid-Loaded Inhomogeneous Elastic Shells of Arbitrary Shape.  _Journal of Vibration and Acoustics_ 115, 384–390 (1993).  https://doi.org/10.1115/1.2930361.
  * [17] Norris, A. N. & Rebinsky, D. A.  Membrane and Flexural Waves on Thin Shells.  _Journal of Vibration and Acoustics_ 116, 457–467 (1994).  https://doi.org/10.1115/1.2930449.
  * [18] Pierce, A. D.  Wave propagation on thin-walled elastic cylindrical shells.  In McCarthy, M. & Hayes, M. (eds.) _Elastic Wave Propagation_ , vol. 35 of _North-Holland Series in Applied Mathematics and Mechanics_ , 205–210 (North-Holland, 1989).  https://doi.org/10.1016/B978-0-444-87272-2.50034-8.
  * [19] Bolotin, V. V.  _Random vibrations of elastic systems_ , vol. 8 (Springer Science & Business Media, 2013).  https://doi.org/10.1007/978-94-017-2842-3.
  * [20] Zhang, X., Liu, G. & Lam, K.  Vibration analysis of thin cylindrical shells using wave propagation approach.  _Journal of Sound and Vibration_ 239, 397–403 (2001).  https://doi.org/10.1006/jsvi.2000.3139.
  * [21] Pellicano, F.  Vibrations of circular cylindrical shells: Theory and experiments.  _Journal of Sound and Vibration_ 303, 154–170 (2007).  https://doi.org/10.1016/j.jsv.2007.01.022.
  * [22] Morsbøl, J., Sorokin, S. & Peake, N.  A WKB approximation of elastic waves travelling on a shell of revolution.  _Journal of Sound and Vibration_ 375, 162–186 (2016).  https://doi.org/10.1016/j.jsv.2016.04.001.
  * [23] Norris, A. N.  Rays, beams and quasimodes on thin shell structures.  _Wave Motion_ 21, 127–147 (1995).  https://doi.org/10.1016/0165-2125(94)00048-A.
  * [24] Green, G. & Fleischmann, R.  Branched flow and caustics in nonlinear waves.  _New Journal of Physics_ 21, 083020 (2019).  https://doi.org/10.1088/1367-2630/ab319b.
  * [25] Degueldre, H., Metzger, J. J., Schultheis, E. & Fleischmann, R.  Channeling of branched flow in weakly scattering anisotropic media.  _Physical Review Letters_ 118, 024301 (2017).  https://doi.org/10.1103/PhysRevLett.118.024301.
  * [26] Degueldre, H.-P.  _Random Focusing of Tsunami Waves_.  Ph.D. thesis, Georg-August-Universitat Gottingen (2015).  http://dx.doi.org/10.53846/goediss-5365.
  * [27] Metzger, J. J.  _Branched Flow and Caustics in Two-Dimensional Random Potentials and Magnetic Fields_.  Ph.D. thesis, Georg-August-University Goettingen (2010).  http://dx.doi.org/10.53846/goediss-2917.

## Acknowledgements

Thanks are due to the University of Southampton for access to Iridis Compute
Cluster; EU’s Horizon 2020 programme under Marie Skłodowska-Curie scheme
(Grant Agreement No. 765636) for financial support; Claus Ibsen, Vestas
aircoil for providing a practical context.

## Author contributions statement

K.J. conceived of and performed the research. K.J. & A.B. analyzed the data.
K.J., N.F. & A.B. discussed the results and wrote the manuscript. A.B. & N.F.
secured the funding and supervised this work.

## Competing interests

The authors declare no competing interests.