transform of the power spectrum, which describes the amplitude of density
fluctuations across distance scales, the 2ptCF is indeed expected to constrain
$\sigma_{8}$ well, even despite higher Poisson noise at the lower density. We
also find that the VPF often results in better constraints on
$\Omega_{\text{M}}$ than the 2ptCF and CiC. The VPF also appears to drive the
bulk of the weak constraints on the SAM SN parameters. Why might this be?

First, let us consider our approach to CiC. CiC in its entirety contains all
information from all orders of correlations (Uhlemann et al., 2020), and might
therefore be expected to give the most robust constraints (e.g. Samushia et
al. 2021). However, due to the limited amount of data the neural network can
take in and still promptly converge on a solution, we cannot use the full
range of our CiC measurements. The input 1D array of clustering should have
fewer than 1000 values, so that our fully-connected neural networks do not
spiral into an unwieldy size with each new layer of neurons. For example, a
fully-connected network of 4 layers with 1,000 neurons with an input of 10,000
elements and output of 5 elements would have:
10,000$\times$1,000${}^{4}\times$5 internal relationships to consider when
optimizing the network. This can be feasible with powerful GPUs and patience,
but the amount of time it takes to train grows very quickly. In our work, even
sampling to cells with fewer than 50 galaxies ($n=50$ in the terminology of
Table LABEL:table:CiCRadii) for 10 distance scales between 1 and 40 cMpc
yields 500 data points for a single redshift.

The default of our current setup uses four redshifts between $0<z<1$, and in
this scenario we instead sample CiC at three distance scales where $n=0-50$,
taking care to choose distance scales that show the differences between the
simulations, as well as aiming to sample evenly across approximately $5-40$
cMpc. Additionally, taking distance scales that are far apart helps reduce the
correlation between measurements (e.g. see Gangolli et al. 2021 for an
examination of correlations within VPF measurements, and Uhlemann et al. 2020
for CiC). Table LABEL:table:CiCRadii fully describes how we sampled a small
part of the full CiC distribution due to computational limitations, both for
our default four-redshift networks, as well as when focusing on one redshift
at a time in §4.3 and 5.2. Though we have sampled a small part of the full CiC
distribution due to computational requirements, we still find the CiC inputs
offer competitive constraints and sometimes outperform the 2ptCF.161616Future
work with the clustering analyzed here could, for example, leverage
convolutional neural networks in order to use the full count-in-cells
distributions for parameter estimation: e.g. as done with CAMELS CMD in
Villaescusa-Navarro et al. (2021c).

We note here the relevant study of Wang et al. (2019) that reached similar
conclusions within a different framework. Using the ‘decorated’ HOD (Hearin et
al., 2016) that includes galaxy assembly bias, they used a complete Fisher
matrix analysis to probe how robustly various clustering statistics constrain
assembly bias. They compared the projected 2ptCF (which projects the 3D 2ptCF
into the dimensions of galaxy observations in RA, Dec, and degrees); the VPF;
the galaxy-galaxy lensing signal; and variations on CiC such as count-in-
cylinders and -annuli and probability distributions of them. Their work
strongly motivates including CiC statistics alongside the popular projected
2ptCF and lensing signal for efficient constraints on galaxy assembly bias.
Specifically, the VPF was good at constraining the number of central galaxies
(not surprisingly, since it is a binary statistic that finds only empty test
spheres, and will be less sensitive to satellites), while the varied CiC
refined the number of centrals and satellites well.

Our handling of CiC may give an explanation for why the VPF performs well. The
VPF is the 0${}^{\text{th}}$ moment of counts-in-cells, the $n=0$ measurement.
Because it yields a single value at each distance scale, we are able to
include many more distance scales when preparing our data for a neural
network. Therefore, the VPF serves as a sampling of the CiC and all the
higher-order moments at many distance scales (White, 1979). Future work may
therefore benefit from including the VPF and similar statistics alongside the
popular 2ptCF (e.g. the $k$-NN statistics proposed by Banerjee & Abel 2020).

##  7 CAMELS-SAM Discussion and Implications

This work is a proof of concept, both for the power of the CAMELS-SAM
simulation suite and the use of SAM-generated galaxy catalogs, galaxy
clustering, and neural networks to constrain cosmology and galaxy formation.
In this section, we discuss implications of our analyses.

###  7.1 Comparison with original hydrodynamic CAMELS

Some readers may wonder if the work carried out in this work can be done to
any extent with the original CAMELS suites. Important to the conception of
CAMELS-SAM was the difficulty of applying galaxy selections across the entire
CAMELS hydrodynamic suite. For example, if one tries to select any objects
with non-zero stellar mass, the SIMBA ‘hump’ will often easily find several
thousand galaxies in the (25 $h^{-1}$ cMpc)3 volumes, while 10 or more percent
of the IllustrisTNG ‘hump’ would struggle to get more than a hundred
galaxies171717As explored in Perez et al. (2021), when, where, and how the VPF
clustering should be measured depends on galaxy density and total covered
volume. This logic can be extended to CiC comfortably, and confirms much of
the common-sense logic of the 2ptCF in the literature.. We found there is no
basic galaxy property selection that yields enough objects for acceptable
Poisson noise and a large enough training sets across both the SIMBA and TNG
humps, even if not down-sampling to a single density. We are therefore unable
to robustly explore how well galaxy clustering marginalizes over the
hydrodynamic astrophysics models to measure cosmology.

Among the most stringent CAMELS selection we are able to make across both
hydrodynamic ‘humps’ is a halo mass cut of M${}_{\text{halo}}>2\times
10^{10}h^{-1}$ M⊙. We make this cut and then randomly down-sample to
$\mathcal{N}=0.064h^{3}$ cMpc-3 (1000 galaxies in each volume). We apply this
to both the IllustrisTNG-DM and SIMBA-DM humps, and update our clustering to
account for the smaller volume (and therefore allowed distance
scales)181818The VPF radii are 10 between with R$=0.8-8\ h^{-1}$ cMpc; the
2ptCF radii are 19 bins whose edges are evenly log-spaced between 0.68 and 9.3
$h^{-1}$ cMpc; and we include CiC to $n=50$ at R $=\\{3.2,5.6,7.2\\}$ cMpc. We
combine the two ‘humps’ for training under the MSE loss function criterion,
yielding a total combined suite of approximately 1400/300/300
training/validation/testing simulations.. Using the clustering of these
CAMELS-DM halos, we find rMSE constraints of 0.031 for $\Omega_{\text{M}}$ and
0.081 for $\sigma_{8}$ (10$\%$ atop the fiducial $\Omega_{\text{M}}=0.3$ and
$\sigma_{8}=0.8$).

Cosmic variance in CAMELS—or the variance due to creating cosmological volumes
of just (25 $h^{-1}$ cMpc)3 for each combination of parameters—results in
noisy clustering statistics, and makes the neural network predictions less
accurate. Our choice to create much larger volumes at lower mass resolution
improved the predictive power of our neural networks, both in terms of
decreasing the effect of cosmic variance and expanding the galaxy selection
and scale of clustering measurements we are able to carry out. This exercise
further motivates the creation of larger CAMELS hydrodynamic wings, and the
development of other techniques such as “next generation” SAMs, which are
specifically designed to emulate the results of specific hydrodynamic
simulations.

###  7.2 CAMELS-SAM in context

Next, we discuss CAMELS-SAM in context with other large simulation suites, and
some future work and possible experiments that CAMELS-SAM enables.

The N-body portion of CAMELS-SAM straddles a unique point between the
limitations of computing power, data storage capacity, and useful scientific
application, especially for machine learning. This part of CAMELS-SAM
comprises over 1000 unique simulations of a moderately large volume, with
moderate mass resolution (sufficient to robustly resolve galaxy properties
relevant to upcoming observational samples), with 100 snapshots saved
throughout $20<z<0$, and sampling of very broad cosmological parameter space
in $\Omega_{\text{M}}$ and $\sigma_{8}$. Individually, other simulation suites
may be comparable to or superior to each of these aspects, but they are
combined in such a way to fill a unique role.

There exist N-body simulation suites that are much larger and/or higher
resolution—e.g. BACCO (Angulo et al., 2021), Aemulus (DeRose et al., 2019),
ABACUS Cosmo/Summit (Garrison et al. 2021; Maksimova et al. 2021), Uchuu
(Ishiyama et al., 2021), Dark Quest (Nishimichi et al., 2019)–but that may not
be as well suited for training neural networks or for running semi-analytic
models. For example, comparable simulation suites often contain significantly
lower numbers of realizations (e.g. between several dozen to several hundred),
which risks providing a small training set191919It is worth noting, however,
that these suites were not created for machine learning training, and they
excel at creating e.g. accurate and robust emulators of various phenomena. We
acknowledge our somewhat unfair comparison.. Next, many of these suites cover
a narrower cosmological parameter space, which can risk neural network results
too tightly focused around the priors (Villaescusa-Navarro et al. 2020a; e.g.
Figure 6 of Ntampaka et al. 2020 with AbacusCosmo, who had to restrict their
cosmological space in $\sigma_{8}$ due to biasing at the edges). Additionally,
many of these suites solved the volume-resolution-data storage balance by
saving a small number of snapshots (often 10-50, with as many as 65 or as few
as 5). Though ideal for their specific science goals, this limits the
possibility of running SAMs to generate galaxies, as they require densely
sampled merger tree histories. Finally, we also point out N-body simulation
suites that do not vary cosmological parameters, but which are well-suited to
study other key features of large-scale structure and cosmology (e.g. Indra
from Falck et al. 2021, several hundred large volumes and many snapshots for
excellent statistics; UNIT from Chuang et al. 2019, several hundred large
volumes at excellent resolution for non-linear statistics).

Most comparable to CAMELS-SAM (beyond the hydrodynamic CAMELS ‘humps’) is the
Quijote project (Villaescusa-Navarro et al., 2020b). The Quijote suite is
unique in that it covers an even broader cosmological parameter space than all
of CAMELS: 7,000 unique models over 6 cosmological parameters, including
massive neutrinos. Quijote also has much larger volumes of 1 ($h^{-1}$ Gpc)3
(though at lower resolution), and an astounding 44,000 simulations. However,
its 5 stored snapshots and lower mass resolution make it unsuitable for
applications using merger trees, and therefore from being the backbone for
CAMELS-SAM. Though dark matter only, the vast Quijote simulation suite may be
able to answer questions this work has not or cannot; for example, probing
larger scale clustering, using more clustering statistics, implementing more
sophisticated statistical tools like the Fisher matrix, expanding the
cosmological models probed, etc.

###  7.3 CAMELS-SAM Data Release & Possibilities

Like the previous CAMELS ‘humps’, CAMELS-SAM was created to serve as a data
set to train machine learning tools to measure and analyze cosmology and
uncertain aspects of astrophysics models. With the addition of the Santa Cruz
SAM-generated galaxies atop large N-body only volumes, CAMELS-SAM offers a
completely unique data set for machine learning. As part of the CAMELS Public
Data Release in Villaescusa-Navarro et al.
(2022)202020https://camels.readthedocs.io/<EMAIL_ADDRESS>we
release:

  * •

The halo catalogues from Rockstar.

  * •

The merger trees generated from Consistent trees.

  * •

The galaxy catalogues from the Santa Cruz SAM.

  * •

Documentation at https://camels-sam.readthedocs.io/

Those wishing to directly analyze the raw simulation snapshots of CAMELS-SAM
should reach out to the CAMELS team. The raw data (full N-body snapshots
across redshifts) has been stored on tape and may be retrieved upon request.

We offer a brief list of ideas of the type of work possible with CAMELS-SAM as
inspiration: machine learning studies of cosmological effects on halo
properties and evolution; creating cosmic void catalogs from both the N-body
only and SC-SAM galaxy catalogs and exploring the effect of cosmology and
astrophysics; how other galaxy observables like the SMF evolve with cosmology
and feedback parameters; and expanding the summary statistics used to
constrain cosmology.

Finally, the dark matter-only data products alone offer many possibilities.
CAMELS-SAM was created specifically to run the Santa Cruz SAM, but other SAMs
could also be implemented within our merger trees. Similarly, the simulations
should offer a fascinating playground for HOD analyses.

Within the realm of our work with CAMELS-SAM, future progress includes:
expanding the clustering statistics and methods applied, varying more
parameters in the SC-SAM, exploring more methods for tightening our
cosmological constraints and marginalization over astrophysics, and
implementing or creating more realistic galaxy selections to yield realistic
mock galaxy catalogs.

##  8 Conclusion

We present CAMELS-SAM, a new and bigger ‘hump’ of the CAMELS project. CAMELS-
SAM is composed of more than 1000 unique N-body only simulations of volume
($100\ h^{-1}$ Mpc)3 and N=$640^{3}$ particles each. The N-body simulations
were generated across a broad cosmological parameter space of
$\Omega_{\text{M}}$=[0.1,0.5] and $\sigma_{8}$=[0.6,1.0], with 100 stored
snapshots between $20\leq z\leq 0$. Each of these N-body simulations has
associated Rockstar halo catalogs and ConsistentTrees merger trees. Finally,
each N-body simulation was run through a unique iteration of the Santa Cruz
SAM for galaxy formation, covering a broad range of the parameters controlling
feedback from massive stars, supernovae, and AGN radio jets.

These halo catalogs, merger trees, and galaxy catalogs have been publicly
released for the community to use for a variety of science applications:
https://camels-sam.readthedocs.io (Villaescusa-Navarro et al., 2022). As a
proof-of-concept for the capabilities of this simulation suite, we used galaxy
clustering statistics to constrain cosmology, marginalize over astrophysics,
and probe astrophysical feedback with simple neural networks.

A brief summary of our work and results:

  * •

We measure and analyze the Void Probability Function, counts-in-cells, and
real-space 3D two-point correlation function of halos and SC-SAM galaxies. We
compare clustering across selections by halo mass, stellar mass, star
formation rate, and specific star formation rate.

  * •

We leverage simple 1D:1D neural networks and a likelihood-free inference
method to measure the mean and standard deviation of each parameter’s marginal
posterior. We leverage CAMELS-SAM’s large suite size for a split of
700/150/150 training/validation/testing simulation sets.

  * •

We explore how the accuracy and precision of our parameter inference varies
with different halo/galaxy selections and density down-samplings, and compare
the choice of combining vs. keeping separate the redshift and clustering
statistics used. §4 focuses on the cosmological parameters
{$\Omega_{\text{M}}$, $\sigma_{8}$}, while §5 focuses on the SC-SAM feedback
parameters.

  * •

The tightest constraints on cosmology we find are with the clustering by dark
matter halo mass. Our cosmological constraints for {$\Omega_{\text{M}}$,
$\sigma_{8}$} find fractional errors of {4.7$\%$, 3$\%$} about their fiducial
values of {0.3, 0.8}, respectively. With the clustering of SC-SAM galaxies by
stellar mass, we predict {$\Omega_{\text{M}}$, $\sigma_{8}$} with errors of
{4.7-6.7$\%$, 3-5$\%$}.

  * •

We find that our other selections based on stellar mass yield predictions for
{$\Omega_{\text{M}}$, $\sigma_{8}$} with errors of {8.6-12.3$\%$, 5.4-8$\%$}.
Selecting on instantaneous star formation rate often yields predictions for
{$\Omega_{\text{M}}$, $\sigma_{8}$} with errors of {10-15.6$\%$,
6.4-12.8$\%$}. Selecting by specific star formation rate often yields
predictions for {$\Omega_{\text{M}}$, $\sigma_{8}$} with errors of
{8.3-9.3$\%$, 4.6-6.9$\%$}.

  * •

We tend to find tighter constraints on both cosmological parameters and
astrophysics parameters when we do _not_ randomly down-sample to a fixed
number density. However, some galaxy property selections perform comparably
well after down-sampling.

  * •

We find better constraints when we combine clustering statistics from several
redshifts. When using one redshift at a time, we do not find evidence for
strong redshift dependence on the quality of the constraints.

  * •

In §6, we compare the constraints that each clustering statistics finds
independently. When comparing each of the clustering statistics–VPF, CiC,
2ptCF–for various mass-threshold clustering samples, we find that all the
statistics find similar constraints for $\sigma_{8}$.

  * •

We find that CiC and VPF often slightly out-perform the 2ptCF alone for
$\Omega_{\text{M}}$. We also find the VPF often drives the bulk of the
constraints found on the SC-SAM parameters.

Finally, these are key implications of our work and CAMELS-SAM:

  * •

Our work contributes to the growing literature that advocates the inclusion of
count-in-cells, the VPF, and related statistics alongside the 2ptCF for
constraining cosmology (e.g. Wang et al. 2019; Uhlemann et al. 2020; Samushia
et al. 2021).

  * •

This is the first work to leverage machine learning and adjacent tools to not
only constrain parameters of a SAM (e.g. van Daalen et al. 2016), but to
simultaneously constrain cosmology and improve the information a neural
network is able to learn from galaxy clustering.

  * •

Our work with CAMELS-SAM includes smaller non-linear scales than most have
probed for cosmology (reaching R$>$1.1-1.6 cMpc, or $k_{\text{max}}<5.85-8.15\
h$ Mpc-1).

  * •

By implementing a robust model for galaxies and their complex astrophysics,
our neural networks are able to learn the unique ways that the SC-SAM affects
galaxy bias, and significantly improve constraints on cosmology than if we
were to simply assume galaxies follow the clustering of dark matter halos.

  * •

The halo products from CAMELS-SAM inhabit a unique position in simulation
volume and resolution, excellent snapshot and redshift coverage, number of
simulations, and vast cosmological parameter space, ideal for several
applications across machine learning, galaxy modeling, and dark matter
analysis.

This work was supported by the Flatiron Institute’s Center for Computational
Astrophysics Pre-Doctoral Program for the fall of 2020. CAMELS-SAM was run,
stored, and remains accessible through the Flatiron Institute Scientific
Computing Hub. DAA was supported in part by NSF grants AST-2009687 and
AST-2108944. We gratefully thank M. Hirschmann for valuable conversations and
explorations about simulating galaxies for clustering early in the project; R.
Angulo for helpful discussion about the BACCO suite and its cosmological
rescaling methods; A. Banerjee for fascinating discussion about their k-NN
method; and L. Garrison for support and helpful context around CORRFUNC and
the Abacus code. We thank P. Behroozi for creating and sharing ROCKSTAR and
ConsistentTrees. We also thank the Flatiron Institute’s Scientific Computing
team for helping make this project feasible and enduring, and especially to J.
Tischio for support with the enormous data space.

## References

  * Abbott et al. (2018) Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526, doi: 10.1103/PhysRevD.98.043526
  * Akiba et al. (2019) Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, in Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining
  * Alimi et al. (2012) Alimi, J.-M., Bouillot, V., Rasera, Y., et al. 2012, arXiv e-prints, arXiv:1206.2838.  https://arxiv.org/abs/1206.2838
  * Amendola et al. (2018) Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Reviews in Relativity, 21, 2, doi: 10.1007/s41114-017-0010-3
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar, D., Faucher-Giguère, C.-A., Kereš, D., et al. 2017, MNRAS, 470, 4698, doi: 10.1093/mnras/stx1517
  * Angulo et al. (2021) Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869, doi: 10.1093/mnras/stab2018
  * Aricò et al. (2021) Aricò, G., Angulo, R. E., Contreras, S., et al. 2021, MNRAS, 506, 4070, doi: 10.1093/mnras/stab1911
  * Baldry et al. (2012) Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621, doi: 10.1111/j.1365-2966.2012.20340.x
  * Banerjee & Abel (2020) Banerjee, A., & Abel, T. 2020, arXiv e-prints, arXiv:2007.13342.  https://arxiv.org/abs/2007.13342
  * Barreira et al. (2021) Barreira, A., Lazeyras, T., & Schmidt, F. 2021, J. Cosmology Astropart. Phys, 2021, 029, doi: 10.1088/1475-7516/2021/08/029
  * Barro et al. (2013) Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104, doi: 10.1088/0004-637X/765/2/104
  * Bayer et al. (2021) Bayer, A. E., Villaescusa-Navarro, F., Massara, E., et al. 2021, ApJ, 919, 24, doi: 10.3847/1538-4357/ac0e91
  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18, doi: 10.1088/0004-637X/763/1/18
  * Bernardi et al. (2013) Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697, doi: 10.1093/mnras/stt1607
  * Blanton & Moustakas (2009) Blanton, M. R., & Moustakas, J. 2009, ARA&A, 47, 159, doi: 10.1146/annurev-astro-082708-101734
  * Borrow et al. (2020) Borrow, J., Anglés-Alcázar, D., & Davé, R. 2020, MNRAS, 491, 6102, doi: 10.1093/mnras/stz3428
  * Brammer et al. (2011) Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2011, ApJ, 739, 24, doi: 10.1088/0004-637X/739/1/24
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Calette et al. (2018) Calette, A. R., Avila-Reese, V., Rodríguez-Puebla, A., Hernández-Toledo, H., & Papastergis, E. 2018, Rev. Mexicana Astron. Astrofis., 54, 443.  https://arxiv.org/abs/1803.07692
  * Calzetti (2013) Calzetti, D. 2013, Star Formation Rate Indicators, ed. J. Falcón-Barroso & J. H. Knapen, 419
  * Catinella et al. (2018) Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875, doi: 10.1093/mnras/sty089
  * Chuang et al. (2019) Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48, doi: 10.1093/mnras/stz1233
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conroy et al. (2005) Conroy, C., Coil, A. L., White, M., et al. 2005, ApJ, 635, 990, doi: 10.1086/497682
  * Contreras et al. (2020) Contreras, S., Angulo, R. E., Zennaro, M., Aricò, G., & Pellejero-Ibañez, M. 2020, MNRAS, 499, 4905, doi: 10.1093/mnras/staa3117
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  * Davidzon et al. (2018) Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107, doi: 10.3847/1538-4357/aaa19e
  * DeRose et al. (2019) DeRose, J., Wechsler, R. H., Tinker, J. L., et al. 2019, ApJ, 875, 69, doi: 10.3847/1538-4357/ab1085
  * Duncan et al. (2014) Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960, doi: 10.1093/mnras/stu1622
  * Ellis (2008) Ellis, R. S. 2008, Observations of the High Redshift Universe, ed. A. Loeb, A. Ferrara, & R. S. Ellis, 259–364, doi: 10.1007/978-3-540-74163-3_3
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455, doi: 10.1146/annurev-astro-081811-125521
  * Falck et al. (2021) Falck, B., Wang, J., Jenkins, A., et al. 2021, MNRAS, 506, 2659, doi: 10.1093/mnras/stab1823
  * Förster Schreiber & Wuyts (2020) Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661, doi: 10.1146/annurev-astro-032620-021910
  * Gabrielpillai et al. (2021) Gabrielpillai, A., Somerville, R. S., Genel, S., et al. 2021, arXiv e-prints, arXiv:2111.03077.  https://arxiv.org/abs/2111.03077
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41, doi: 10.1111/j.1365-2966.2005.09321.x
  * Gangolli et al. (2021) Gangolli, N., D’Aloisio, A., Nasir, F., & Zheng, Z. 2021, MNRAS, 501, 5294, doi: 10.1093/mnras/staa3843
  * Garrison et al. (2021) Garrison, L. H., Eisenstein, D. J., Ferrer, D., Maksimova, N. A., & Pinto, P. A. 2021, MNRAS, 508, 575, doi: 10.1093/mnras/stab2482
  * Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175, doi: 10.1093/mnras/stu1654
  * Gonzalez-Perez et al. (2018) Gonzalez-Perez, V., Comparat, J., Norberg, P., et al. 2018, MNRAS, 474, 4024, doi: 10.1093/mnras/stx2807
  * Guo et al. (2013) Guo, Q., White, S., Angulo, R. E., et al. 2013, MNRAS, 428, 1351, doi: 10.1093/mnras/sts115
  * Guo et al. (2011) Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101, doi: 10.1111/j.1365-2966.2010.18114.x
  * Habouzit et al. (2020) Habouzit, M., Pisani, A., Goulding, A., et al. 2020, MNRAS, 493, 899, doi: 10.1093/mnras/staa219
  * Hadzhiyska et al. (2021a) Hadzhiyska, B., Liu, S., Somerville, R. S., et al. 2021a, MNRAS, 508, 698, doi: 10.1093/mnras/stab2564
  * Hadzhiyska et al. (2021b) Hadzhiyska, B., Tacchella, S., Bose, S., & Eisenstein, D. J. 2021b, MNRAS, 502, 3599, doi: 10.1093/mnras/stab243
  * Hahn & Villaescusa-Navarro (2021) Hahn, C., & Villaescusa-Navarro, F. 2021, J. Cosmology Astropart. Phys, 2021, 029, doi: 10.1088/1475-7516/2021/04/029
  * Hamaus et al. (2016) Hamaus, N., Pisani, A., Sutter, P. M., et al. 2016, Phys. Rev. Lett., 117, 091302, doi: 10.1103/PhysRevLett.117.091302
  * Hamilton (2001) Hamilton, A. J. S. 2001, MNRAS, 322, 419, doi: 10.1046/j.1365-8711.2001.04137.x
  * Hassan et al. (2021) Hassan, S., Villaescusa-Navarro, F., Wandelt, B., et al. 2021, arXiv e-prints, arXiv:2110.02983.  https://arxiv.org/abs/2110.02983
  * Hearin et al. (2016) Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016, MNRAS, 460, 2552, doi: 10.1093/mnras/stw840
  * Ishiyama et al. (2021) Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210, doi: 10.1093/mnras/stab1755
  * Jeffrey et al. (2022) Jeffrey, N., Boulanger, F., Wandelt, B. D., et al. 2022, MNRAS, 510, L1, doi: 10.1093/mnrasl/slab120
  * Jeffrey & Wandelt (2020) Jeffrey, N., & Wandelt, B. D. 2020, arXiv e-prints, arXiv:2011.05991.  https://arxiv.org/abs/2011.05991
  * Jiang et al. (2021) Jiang, F., Dekel, A., Freundlich, J., et al. 2021, MNRAS, 502, 621, doi: 10.1093/mnras/staa4034
  * Khostovan et al. (2018) Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2018, MNRAS, 478, 2999, doi: 10.1093/mnras/sty925
  * Khostovan et al. (2020) Khostovan, A. A., Malhotra, S., Rhoads, J. E., et al. 2020, MNRAS, 493, 3966, doi: 10.1093/mnras/staa175
  * Kirby et al. (2011) Kirby, E. N., Lanfranchi, G. A., Simon, J. D., Cohen, J. G., & Guhathakurta, P. 2011, ApJ, 727, 78, doi: 10.1088/0004-637X/727/2/78
  * Kokron et al. (2021) Kokron, N., DeRose, J., Chen, S.-F., White, M., & Wechsler, R. H. 2021, MNRAS, 505, 1422, doi: 10.1093/mnras/stab1358
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Landy & Szalay (1993) Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64, doi: 10.1086/172900
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473, doi: 10.1086/309179
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Maksimova et al. (2021) Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, MNRAS, 508, 4017, doi: 10.1093/mnras/stab2484
  * Massara et al. (2021) Massara, E., Villaescusa-Navarro, F., Ho, S., Dalal, N., & Spergel, D. N. 2021, Phys. Rev. Lett., 126, 011301, doi: 10.1103/PhysRevLett.126.011301
  * McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184, doi: 10.1088/0004-637X/764/2/184
  * McCullagh et al. (2017) McCullagh, N., Norberg, P., Cole, S., et al. 2017, arXiv e-prints, arXiv:1705.01988.  https://arxiv.org/abs/1705.01988
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Mead et al. (2021) Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401, doi: 10.1093/mnras/stab082
  * Mobasher et al. (2015) Mobasher, B., Dahlen, T., Ferguson, H. C., et al. 2015, ApJ, 808, 101, doi: 10.1088/0004-637X/808/1/101
  * Moser et al. (2022) Moser, E., Battaglia, N., Nagai, D., et al. 2022, arXiv e-prints, arXiv:2201.02708.  https://arxiv.org/abs/2201.02708
  * Moustakas et al. (2013) Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50, doi: 10.1088/0004-637X/767/1/50
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Naab & Ostriker (2017) Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59, doi: 10.1146/annurev-astro-081913-040019
  * Netzer (2015) Netzer, H. 2015, ARA&A, 53, 365, doi: 10.1146/annurev-astro-082214-122302
  * Nicola et al. (2022) Nicola, A., Villaescusa-Navarro, F., Spergel, D. N., et al. 2022, arXiv e-prints, arXiv:2201.04142.  https://arxiv.org/abs/2201.04142
  * Nishimichi et al. (2019) Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29, doi: 10.3847/1538-4357/ab3719
  * Ntampaka et al. (2020) Ntampaka, M., Eisenstein, D. J., Yuan, S., & Garrison, L. H. 2020, ApJ, 889, 151, doi: 10.3847/1538-4357/ab5f5e
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, in Advances in Neural Information Processing Systems 32, ed. H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, & R. Garnett (Curran Associates, Inc.), 8024–8035
  * Peebles (1980) Peebles, P. J. E. 1980, The large-scale structure of the universe
  * Perez et al. (2021) Perez, L. A., Malhotra, S., Rhoads, J. E., & Tilvi, V. 2021, ApJ, 906, 58, doi: 10.3847/1538-4357/abc88b
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Pisani et al. (2015) Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, Phys. Rev. D, 92, 083531, doi: 10.1103/PhysRevD.92.083531
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Porter et al. (2014) Porter, L. A., Somerville, R. S., Primack, J. R., & Johansson, P. H. 2014, MNRAS, 444, 942, doi: 10.1093/mnras/stu1434
  * Raichoor et al. (2017) Raichoor, A., Comparat, J., Delubac, T., et al. 2017, MNRAS, 471, 3955, doi: 10.1093/mnras/stx1790
  * Repp & Szapudi (2020) Repp, A., & Szapudi, I. 2020, MNRAS, 498, L125, doi: 10.1093/mnrasl/slaa139
  * Riccio et al. (2021) Riccio, G., Małek, K., Nanni, A., et al. 2021, A&A, 653, A107, doi: 10.1051/0004-6361/202140854
  * Rodríguez-Puebla et al. (2017) Rodríguez-Puebla, A., Primack, J. R., Avila-Reese, V., & Faber, S. M. 2017, MNRAS, 470, 651, doi: 10.1093/mnras/stx1172
  * Salvador et al. (2019) Salvador, A. I., Sánchez, F. J., Pagul, A., et al. 2019, MNRAS, 482, 1435, doi: 10.1093/mnras/sty2802
  * Samushia et al. (2021) Samushia, L., Slepian, Z., & Villaescusa-Navarro, F. 2021, MNRAS, 505, 628, doi: 10.1093/mnras/stab1199
  * Sinha & Garrison (2020) Sinha, M., & Garrison, L. H. 2020, MNRAS, 491, 3022, doi: 10.1093/mnras/stz3157
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Somerville et al. (2008) Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481, doi: 10.1111/j.1365-2966.2008.13805.x
  * Somerville et al. (2015) Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337, doi: 10.1093/mnras/stv1877
  * Somerville & Primack (1999) Somerville, R. S., & Primack, J. R. 1999, MNRAS, 310, 1087, doi: 10.1046/j.1365-8711.1999.03032.x
  * Somerville et al. (2021) Somerville, R. S., Olsen, C., Yung, L. Y. A., et al. 2021, MNRAS, 502, 4858, doi: 10.1093/mnras/stab231
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Steinhardt & Speagle (2014) Steinhardt, C. L., & Speagle, J. S. 2014, ApJ, 796, 25, doi: 10.1088/0004-637X/796/1/25
  * Sugiyama et al. (2020) Sugiyama, S., Takada, M., Kobayashi, Y., et al. 2020, Phys. Rev. D, 102, 083520, doi: 10.1103/PhysRevD.102.083520
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253, doi: 10.1086/191823
  * Szewciw et al. (2021) Szewciw, A. O., Beltz-Mohrmann, G. D., Berlind, A. A., & Sinha, M. 2021, arXiv e-prints, arXiv:2110.03701.  https://arxiv.org/abs/2110.03701
  * Uhlemann et al. (2020) Uhlemann, C., Friedrich, O., Villaescusa-Navarro, F., Banerjee, A., & Codis, S. 2020, MNRAS, 495, 4006, doi: 10.1093/mnras/staa1155
  * Vakili et al. (2020) Vakili, M., Hoekstra, H., Bilicki, M., et al. 2020, arXiv e-prints, arXiv:2008.13154.  https://arxiv.org/abs/2008.13154
  * van Daalen et al. (2016) van Daalen, M. P., Henriques, B. M. B., Angulo, R. E., & White, S. D. M. 2016, MNRAS, 458, 934, doi: 10.1093/mnras/stw405
  * Villaescusa-Navarro et al. (2020a) Villaescusa-Navarro, F., Wandelt, B. D., Anglés-Alcázar, D., et al. 2020a, arXiv e-prints, arXiv:2011.05992.  https://arxiv.org/abs/2011.05992
  * Villaescusa-Navarro et al. (2020b) Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020b, ApJS, 250, 2, doi: 10.3847/1538-4365/ab9d82
  * Villaescusa-Navarro et al. (2021a) Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021a, ApJ, 915, 71, doi: 10.3847/1538-4357/abf7ba
  * Villaescusa-Navarro et al. (2021b) Villaescusa-Navarro, F., Genel, S., Angles-Alcazar, D., et al. 2021b, arXiv e-prints, arXiv:2109.10360.  https://arxiv.org/abs/2109.10360
  * Villaescusa-Navarro et al. (2021c) —. 2021c, arXiv e-prints, arXiv:2109.10915.  https://arxiv.org/abs/2109.10915
  * Villaescusa-Navarro et al. (2022) Villaescusa-Navarro, F., Genel, S., Anglés-Alcázar, D., et al. 2022, arXiv e-prints, arXiv:2201.01300.  https://arxiv.org/abs/2201.01300
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031, doi: 10.1093/mnras/stt1789
  * Walcher et al. (2011) Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1, doi: 10.1007/s10509-010-0458-z
  * Walsh & Tinker (2019) Walsh, K., & Tinker, J. 2019, MNRAS, 488, 470, doi: 10.1093/mnras/stz1351
  * Wang et al. (2019) Wang, K., Mao, Y.-Y., Zentner, A. R., et al. 2019, MNRAS, 488, 3541, doi: 10.1093/mnras/stz1733
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435, doi: 10.1146/annurev-astro-081817-051756
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Wen et al. (2020) Wen, D., Kemball, A. J., & Saslaw, W. C. 2020, ApJ, 890, 160, doi: 10.3847/1538-4357/ab6d6f
  * White (1979) White, S. D. M. 1979, MNRAS, 186, 145, doi: 10.1093/mnras/186.2.145
  * Xu et al. (2021) Xu, X., Kumar, S., Zehavi, I., & Contreras, S. 2021, MNRAS, 507, 4879, doi: 10.1093/mnras/stab2464
  * Yung et al. (2019a) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019a, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241
  * Yung et al. (2020) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2020, MNRAS, 496, 4574, doi: 10.1093/mnras/staa1800
  * Yung et al. (2019b) Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, MNRAS, 490, 2855, doi: 10.1093/mnras/stz2755
  * Zehavi et al. (2005) Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2005, ApJ, 630, 1, doi: 10.1086/431891
  * Zhai et al. (2019a) Zhai, Z., Benson, A., Wang, Y., Yepes, G., & Chuang, C.-H. 2019a, MNRAS, 490, 3667, doi: 10.1093/mnras/stz2844
  * Zhai et al. (2019b) Zhai, Z., Tinker, J. L., Becker, M. R., et al. 2019b, ApJ, 874, 95, doi: 10.3847/1538-4357/ab0d7b
  * Zhang et al. (2020) Zhang, G., Li, Z., Liu, J., et al. 2020, Phys. Rev. D, 102, 083537, doi: 10.1103/PhysRevD.102.083537
  * Zou et al. (2019) Zou, H., Gao, J., Zhou, X., & Kong, X. 2019, ApJS, 242, 8, doi: 10.3847/1538-4365/ab1847

##  Appendix A Verifying our SC-SAM runs against observations and
IllustrisTNG300

Figure 12: Verification of our fiducial model and set-up of the SC-SAM at
$z=0$. The 5 CAMELS-SAM ‘cosmic variance’ (CV) galaxy catalogs (red dotted
lines) are compared against relevant $z\sim 0$ observations (light green
shapes) and the default SC-SAM run through IllustrisTNG300-1 (dark gray solid
line; N=25003) and IllustrisTNG300-2 (light gray solid line; N=12503) dark
matter-only volumes. Top row: stellar mass function (left), stellar mass-halo
mass relationship (right). Bottom row: cold gas fraction vs. stellar mass for
disk-dominated galaxies (left), stellar metallicity vs. stellar mass (middle),
black hole mass vs. bulge mass (right). The observational data (green) used
for the calibration of the SC-SAM are detailed in Somerville et al. (2015) and
Yung et al. (2019a). The CV simulations were each generated with a different
random initial seed, and all other parameters set to the fiducial value of
{$\Omega_{\text{M}}$, $\sigma_{8}$, A${}_{\text{SN1}}$, A${}_{\text{SN2}}$,
A${}_{\text{AGN}}$}={0.3, 0.8, 1, 0, 1,}.

To confirm that our implementation of the SC-SAM is running correctly, we
computed various useful summary statistics for the CV simulations and compared
them to larger-volume SC-SAM runs and observations. In Figure 12, we compare
the SAM output of the CV simulations, the SAM output from the two highest
resolution IllustrisTNG300 volumes, and various local universe observations
for: the stellar mass function, the stellar mass-halo mass relation, the cold
gas fraction vs. stellar mass, stellar metallicity-stellar mass, and black
hole mass-bulge mass relationships. The local universe observations were
compiled in Yung et al. (2019a) and include: {Bernardi et al. 2013; Moustakas
et al. 2013; Baldry et al. 2012} for the stellar mass function; {Rodríguez-
Puebla et al. 2017} for stellar mass-halo mass; {Catinella et al. 2018;
Calette et al. 2018} for cold gas-stellar mass for disk-dominated galaxies
(bulge-to-stellar mass ratio is less than 40%); Gallazzi et al. 2005; Kirby et
al. 2011 for stellar metallicity vs. stellar mass; and McConnell & Ma 2013;
Kormendy & Ho 2013 for black hole mass vs. bulge mass relationships.

The great consistency with the larger IllustrisTNG300 SAM catalogs, and all
SAM’s agreement with local observations, support our use of the SC-SAM in this
volume and resolution. We direct readers to the Figure 3 in Gabrielpillai et
al. (2021) for more details for how the SC-SAM was run over IllustrisTNG300-1
(the highest resolution 300 volume), additional comparisons to even higher
resolution volumes, and error bounds for most of the relations (removed here
to minimize visual clutter). We note that IllustrisTNG used the cosmology of
Planck Collaboration et al. (2016, Table 4), with $\Omega_{M}=0.3089$,
$\sigma_{8}=0.8159$, and $h=0.6774$. Though slightly different from the
‘fiducial’ cosmology we used for the CV simulations, the SC-SAM still finds
excellent agreement for these galaxy observables without retuning.

Next, we examine the behavior of the “1P” suite for the SC-SAM parameters
across these galaxy relationships. In CAMELS-SAM, the “1P” set is 12 galaxy
catalogs that probe the minimum and maximum values of the SC-SAM parameters in
our suite. We ran the SC-SAM with one of the SAM pre-factors at a time set to
the minimum or maximum value atop the ‘fiducial’ cosmology simulations
CV$\\_$0 and CV$\\_$1\. This creates three pairs of SC-SAM galaxy catalogs for
each of the 2 N-body simulations. Figures 13, 14, and 15 show key SC-SAM
verification relationships for the 1P galaxy catalogs for A${}_{\text{SN1}}$,
A${}_{\text{SN2}}$, and A${}_{\text{AGN}}$ respectively.

The stellar mass function and stellar mass-halo mass relation shows great
sensitivity to A${}_{\text{SN1}}$ and A${}_{\text{SN2}}$, each with unique
effects. Lower values of A${}_{\text{SN1}}$ (weaker feedback) lift and
increase the low-mass cutoff of the SMF slightly, while lower
A${}_{\text{SN2}}$ much more strongly sharpens the shape and increases the
mass of most galaxies. High values of A${}_{\text{SN1}}$ very strongly lowers
the normalization of the stellar mass-halo mass relationship everywhere,
whereas high A${}_{\text{SN2}}$ makes the low-mass end of the stellar mass-
halo mass relationship steeper. Finally, A${}_{\text{SN1}}$ mostly shifts the
stellar metallicity-stellar mass relationship up (weak feedback) and down
(strong feedback), while A${}_{\text{SN2}}$ affects the slope at all but the
highest masses. Both A${}_{\text{SN}}$ parameters affect the low-mass half of
the cold gas fraction vs. stellar mass relationship in similar ways: lower
values create a sharper drop and lower valleys.

A${}_{\text{AGN}}$, on the other hand, has much more subtle effects. It has
nearly no effect on stellar metallicity vs. stellar mass relationship, and
shows only very mild effects on the SMF and black hole vs. bulge mass
relationship at the highest stellar masses. Unlike the A${}_{\text{SN}}$
parameters, its effect on the relationship of cold gas fraction vs. stellar
mass is on the higher mass end. It does, however, show strong effects on the
stellar mass-halo mass relationship on the right/higher-mass half of the
‘mountain’, suppressing the stellar mass values as its effect is strengthened.

Figure 13: Key galaxy relationships for the 1P pair of
A${}_{\text{SN1}}=0.25\times\epsilon_{\text{SN0}}$ (cyan) and
A${}_{\text{SN1}}=4.0\times\epsilon_{\text{SN0}}$ (magenta). See Fig. 12 for
explanation of plotted quantities.  Figure 14: Key galaxy relationships for
the 1P pair of A${}_{\text{SN2}}=-2+\alpha_{\text{rh}}$ (cyan) and
A${}_{\text{SN2}}=+2+\alpha_{\text{rh}}$ (magenta). See Fig. 12 for
explanation of plotted quantities. Figure 15: Key galaxy relationships for the
1P pair of A${}_{\text{AGN}}=0.25\times\kappa_{\text{radio}}$ (cyan) and
A${}_{\text{AGN}}=4.0\times\kappa_{\text{radio}}$ (magenta). See Fig. 12 for
explanation of plotted quantities.

##  Appendix B Tests for neural network accuracy, continued

As discussed in §3.2.3, part of our verification of our neural network’s
accuracy under the likelihood-free inference involves examining the
distribution of $\chi_{j}(\theta_{i})$ for the test set. This confirmation of
the LFI is nearly identical to that in Figure 3 of Jeffrey et al. (2022). The
authors developed the moment networks framework at the heart of our LFI loss
in Jeffrey & Wandelt (2020), and apply it in the context of detecting the
primordial B-mode of the Cosmic Microwave Background. Even with a model with
more than $10^{5}$ parameters, the moment networks gave good results with
accurately estimated errors. We apply the same test to probe how
underestimated or inaccurate our errors may be.

Figure 16 shows the distribution of $\chi_{j}$ (for each $j$ simulation in the
test set) for two of our clustering selections’ neural network predictions
each parameter $\theta_{i}$. We show the top 5 best-performing neural networks
as a unique color. As discussed in §3.2, the top few best-performing neural
networks are all quite similar in overall behavior, so this highlights the
general behavior of the LFI predictions for the given selections. When the LFI
loss leads to accurate predictions, the distribution of $\chi(\theta_{i})$
will approximate a Gaussian. Unconstrained parameters will appear flat within
a similar range or peak at one of the extreme ends.

Figure 16a shows the $\chi_{j}$ distribution for the marginal posteriors given
‘all’ clustering of dark matter halos with mass greater than $2\times 10^{11}$
M⊙ and randomly down-sampled to 0.005 $h^{3}$ cMpc-3. The cosmological
constraints tend to be very good, so their $\chi_{j}$ distributions are
consistent with a Gaussian distribution centered at zero and with variance of
one. The constraints on the SC-SAM astrophysical parameters are nonexistent
(just around the mean of the prior and with large errors), which show up as
$\chi_{j}$ distributions that are flat along the range or that peak at the
Gaussian tails. Figure 16b shows the $\chi_{j}$ distribution for the marginal
posteriors given ‘all’ clustering of SC-SAM galaxies with stellar mass greater
than $1\times 10^{9}$ M⊙ and randomly down-sampled to 0.005 $h^{3}$ cMpc-3.
This selection gives good constraints on the cosmological parameters and
A${}_{\text{SN2}}$. Even though the constraints are not as good on
A${}_{\text{SN2}}$ (i.e. spread around the 1:1 relationship with large error
bars), the errors are not underestimated.

(a) ‘All’ clustering of dark matter halos with mass greater than $2\times
10^{11}$ M⊙ and randomly down-sampled to 0.005 $h^{3}$ cMpc-3.

(b) ‘All’ clustering of SC-SAM galaxies with stellar mass greater than
$1\times 10^{9}$ M⊙ and randomly down-sampled to 0.005 $h^{3}$ cMpc-3.

Figure 16: The $\chi_{i}$ distribution across the test set for the top 5 best
performing networks (each in a unique color). The dotted grey curve is a
Gaussian whose center is at zero (black vertical line) and has a variance of
one.

##  Appendix C Additional Figures: Testing Effect of Redshift on Constraints

Throughout this work, our default input to our neural networks has been
clustering statistics from $z=\\{0.0,0.1,0.5,1.0\\}$ combined together. This,
for example, mimics possible future experiments leveraging similarly selected
galaxy populations at different redshifts. However, in §4.3 and 5.2, we
considered how might our constraints changed if made only at a single
redshift?

We describe our clustering methodology in §3.1, and specifically note Table
LABEL:table:CiCRadii for the slight adjustments we made to what clustering we
give to the single redshift neural network. We note that a more fair
comparison would have followed the example of the clustering statistic tests
in §6, where the only adjustment we made was splitting the data up by
statistic. However, we expand the CiC distributions given to the neural
network at each redshift to answer: how good might our constraints be if we
focus in on a single redshift, and give a neural network as much data as it
can handle?

We share two representative examples of these neural networks in this
appendix, both down-sampled to a density of 0.005 $h^{3}$ cMpc-3: cosmology-
only constraints for a halo mass selection of $2\times 10^{11}$ M⊙ (Figures 17
and 18), and constraints on all 5 parameters at once for a stellar mass
selection of $1\times 10^{9}$ M⊙ (Figure 19).

(a) $z=1.0$.

(b) $z=0.5$.

(c) $z=0.1$.

(d) $z=0.0$.

Figure 17: Examining redshift dependence on constraints for the the
cosmological parameter $\Omega_{\text{M}}$, using the clustering of halos with
mass greater than $2\times 10^{11}$ M⊙ down-sampled to a density 0.005 $h^{3}$
cMpc-3. See Table LABEL:table:CiCRadii for details for what exact distance
scales and values were used for CiC for ‘all’ clustering measured for this
training.

(a) $z=1.0$

(b) $z=0.5$

(c) $z=0.1$

(d) $z=0.0$

Figure 18: Similar to Figure 17, but for $\sigma_{8}$.

(a) Constraints with $z=1.0$ clustering.

(b) Constraints with $z=0.5$ clustering.

(c) Constraints with $z=0.1$ clustering.

(d) Constraints with $z=0$ clustering.

Figure 19: Examining redshift dependence on cosmological constraints for the
clustering statistics of SAM galaxies with stellar mass greater than $2\times
10^{10}$ M⊙, down-sampled to a density of 0.001 $h^{3}$ cMpc-3. See Table
LABEL:table:CiCRadii for details for what exact distance scales and values
were used for CiC for ‘all’ clustering measured for this training.

##  Appendix D Additional Figures: Comparing Constraints between Clustering
Statistics

In §6, we compare the constraints that each independent clustering statistic
is able to find for our five parameters. Though we tested several selections,
we include two representative examples in this Appendix: cosmology-only
constraints on $\Omega_{\text{M}}$ and $\sigma_{8}$ for a halo mass selection
of $2\times 10^{11}$ M⊙ down-sampled to a density of 0.005 $h^{3}$ cMpc-3
(Figures 20 and 21), and constraints on all 5 parameters at once for a stellar
mass selection of $1\times 10^{9}$ M⊙ (Figure 22). Outliers, marked in red
stars in Figure 22, are simulations in the test set whose ‘$Z$-value’ (Eq. 12)
are greater than 6, and are excluded when calculating the rMSE (see §3.2.3 for
more details).

Finally, we remind readers that the true parameter distributions in the test
set only appear slightly skewed to lower values for the A${}_{\text{SN1}}$ and
A${}_{\text{AGN}}$ parameters, due to their original generation in logarithmic
space. This effect would disappear if plotted in log-scale, though we choose
to keep all scales linear for consistency with plots in Villaescusa-Navarro et
al. (2021a).

(a) 2-point correlation function.

(b) Count-in-Cells.

(c) Void Probability Function.

(d) ‘All’ clustering statistics.

Figure 20: Which clustering statistic is best at constraining cosmology
through $\Omega_{\text{M}}$? We compare the constraints on $\Omega_{\text{M}}$
found by different clustering statistics based on the clustering of dark
matter halos with mass greater than $2\times 10^{11}$ M⊙, randomly sampled to
a density of 0.005 $h^{3}$ cMpc-3. We combine the clustering at
$z=\\{0.0,0.1,0.5,1.0\\}$. The 2ptCF (a) is measured between $1.1<R<36.1$
cMpc; CiC (b) is measured at at $R=16.0,22.4,28.8$ cMpc; and the VPF (c) is
measured between $1.6<R<40$ cMpc. We combine ‘all’ these statistics (d) for
the best constraints. Detailed quantitative comparisons can be found in Table
LABEL:table:ClusteringStats_CosmoConstraints.

(a) 2-point correlation function.

(b) Count-in-Cells.

(c) Void Probability Function.

(d) ‘All’ clustering statistics.

Figure 21: Which clustering statistic is best at constraining $\sigma_{8}$?
The same clustering and set-up as in Figure 20, instead for $\sigma_{8}$.

(a) When training a neural network on the two-point correlation function
between $1.1<R<36.1$ cMpc at $z=\\{0.0,0.1,0.5,1.0\\}$.

(b) When training a neural network on the count-in-cells distribution at
$R=16.0,22.4,28.8$ cMpc at $z=\\{0.0,0.1,0.5,1.0\\}$.

(c) When training a neural network on the Void Probability Function between
$1.6<R<40$ cMpc at $z=\\{0.0,0.1,0.5,1.0\\}$.

(d) When training a neural network on ‘all’ clustering statistics described
above together.

Figure 22:  Which clustering statistic is best at constraining cosmology and
astrophysical parameters? We compare the constraints on $\Omega_{\text{M}}$,
$\sigma_{8}$, A${}_{\text{SN1}}$, A${}_{\text{SN2}}$, and A${}_{\text{AGN}}$
found by each of the clustering statistics that we use: (a) the 2ptCF , (b)
CiC , (c) the VPF, and (d) all combined. We use the respective clustering of
SAM galaxies with stellar mass greater than $1\times 10^{9}$ M⊙, down-sampled
to a density of 0.005 $h^{3}$ cMpc-3. Detailed quantitative comparisons can be
found in Tables LABEL:table:ClusteringStats_CosmoConstraints and
LABEL:table:ClusteringStats_AsnConstraints.


# First-order transitions in spin chains coupled to quantum baths

C. A. Perroni1, A. De Candia1,2, V. Cataudella1, R. Fazio1,3,4, and G. De
Filippis1 1 CNR-SPIN and Physics Department E. Pancini - Università degli
Studi di Napoli Federico II,
Complesso Universitario Monte S. Angelo - Via Cintia - I-80126 Napoli, Italy
2 INFN, Sezione di Napoli - Complesso Universitario di Monte S. Angelo -
I-80126 Napoli, Italy
3 ICTP, Strada Costiera 11, I-34151 Trieste, Italy
4 NEST, Istituto Nanoscienze-CNR, I-56126 Pisa, Italy

###### Abstract

We show that tailoring the dissipative environment allows to change the
features of continuous quantum phase transitions and, even, induce first order
transitions in ferromagnetic spin chains. In particular, using a numerically
exact quantum Monte Carlo method for the paradigmatic Ising chain of one-half
spins in a transverse magnetic field, we find that spin couplings to local
quantum boson baths in the Ohmic regime can drive the transition from the
second to the first order even for a low dissipation strength. Moreover, using
a variational mean-field approach for the treatment of spin-spin and spin-
boson interactions, we point out that phase discontinuities are ascribable to
a dissipation induced effective magnetic field which is intrinsically related
to the bath quantum fluctuations and vanishes for classical baths. The
effective field is able to switch the sign of the magnetization along the
direction of spin-spin interactions. The results can be potentially tested in
recent quantum simulators and are relevant for quantum sensing since the spin
system could not only detect the properties of non-classical baths, but also
the effects of weak magnetic fields.

Introduction. A quantum phase transition (QPT) occurs at zero temperature when
quantum fluctuations, tuned through a physical parameter, such as a magnetic
field, are able to induce a sudden change in the ground state of a system with
many-body interactions [1]. The behavior of static quantities at QPT is
typically obtained at the thermodynamic equilibrium by exploiting the quantum-
to-classical mapping, that is the mapping to a classical statistical model
with an additional imaginary-time dimension [1]. Most studied QPT are
continuous [1, 2]. A paradigmatic example is provided by the Ising chain of
one-half spins where the increase of a transverse magnetic field induces the
change from an ordered ferromagnetic to a paramagnetic disordered phase [3].
Attention has recently focused also on first order QPT [2] whose large
sensitivity to external perturbations can be exploited for sensing
applications.

Any quantum system is inevitably coupled to the environment whose interactions
can significantly change physical features. Then the fundamental question is:
how the coupling with its surroundings affects a system close to a QPT? To
this aim two different routes have been proposed in the literature. In the
former one, the steady state reached by the system coupled to Markovian baths
has been investigated starting from a Lindblad master equation [4, 5]. In the
latter one, dissipation is explictly introduced by modelling the environment
as infinite set of harmonic quantum oscillators [5, 6, 7]. Within this
approach, the entire universe (system+environment) is considered at the
thermodynamical equilibrium. Of course, this does not imply any assumption for
the system stationary state. As a further consequence, non-Markovian effects
can be fully included [8]. In the following we will focus our attention on
this second proposal.

The recent realization of programmable spin models with tunable interactions
[9, 10, 11] has stimulated an intense theoretical study of the transverse
field Ising model in the presence of dissipation [12]. In some studies of this
model, the environment is accurately modelled as an infinite set of local
boson baths. Indeed, each spin is coupled to an infinite number of oscillators
giving rise to the well known spin-boson model [6, 7, 13]. In the Ohmic case,
quantum Monte Carlo (QMC) studies have shown that these dissipative mechanisms
can drive QPT even in the case of a single spin system [13, 14, 15]. For many-
body systems, non perturbative properties depend in a crucial way on the
specific coupling with the bath. Indeed, for the quantum Ising chain coupled
to oscillator baths [16, 17, 18], when spin-bath couplings are along the
direction of spin-spin interactions, the dynamical critical exponent changes
while the transition remains continuous.

In this letter, the spin-bath interaction is further tailored with the aim to
induce first order transitions. In particular, we study zero temperature
stationary properties of the transverse field Ising chain with one-half spins
coupled to bosonic local baths in the Ohmic regime through a term which, for a
single spin, is a combination of Jaynes-Cummings and anti-Jaynes-Cummings
interactions [19, 20]. We use the numerically exact QMC method up to the
thermodynamic limit showing that the spin couplings considered in this work
are able to drive the phase transition from the second to the first order also
for a low dissipation strength. Moreover, we develop a semi-analytical zero
temperature variational mean-field (VMF) approach which is in very good
agreement with QMC method clarifying that phase discontinuities are present
only for quantum baths. In fact, the VMF approach highlights the role played
by a dissipation induced effective magnetic field which switches its sign
along the direction of spin-spin interactions with varying the properties of
the bath. Finally, we point out that the strong sensitivity of the spin system
to the properties of non-classical baths at minimal coupling and to the
strength of the magnetic field could be exploited in phenomena relevant for
the realization of small-scale quantum sensors.

The Model. The Hamiltonian is:

$H=H_{S}+H_{B}+H_{SB},$ (1)

where $H_{S}$ describes the spin couplings within the Ising chain in a
transverse field [3]

$H_{S}=\frac{\Delta}{2}\sum_{i=1}^{L}\sigma^{x}_{i}-\frac{J}{4}\sum_{i=1}^{L}\sigma^{z}_{i}\sigma^{z}_{i+1},$
(2)

with the energy $\Delta$ providing the strength of the transverse field along
x-direction, the energy $J$ the exchange couplings of spins along z-direction,
$i=1,..,L$ indicating the $L$ sites of the chain, $\sigma^{x}_{i}$,
$\sigma^{z}_{i}$ the Pauli matrices on each chain site with eigenvalues
$1$,$-1$. In Eq. (1), the Hamiltonian $H_{B}$ describes $L$ local baths, being
each one associated to one of the L sites $i$:

$H_{B}=\sum_{i,k}\hbar\omega_{k}a_{i,k}^{+}a_{i,k},$ (3)

where $a_{i,k}^{+}$ ($a_{i,k}$) creates (destroys) the boson mode $k$ of the
bath at the site $i$. All the local baths are assumed to have the same
frequency spectrum $\omega_{k}$ (independent of the site $i$). Finally, the
spin-bath coupling combines Jaynes-Cummings (rotating) and anti-Jaynes-
Cummings (counter-rotating) interactions [19, 20] through the dimensionless
bath parameter $\gamma$:

$H_{SB}=\sum_{i=1}^{L}\sum_{k}\lambda_{k}\left\\{a_{i,k}\left[\gamma\sigma_{i}^{-}+(1-\gamma)\sigma_{i}^{+}\right]+\text{h.c.}\right\\},$
(4)

where the bath spectral function $F(\omega)$ is defined in terms of the
couplings $\lambda_{k}$:
$F(\omega)=\sum_{k}\lambda_{k}^{2}\delta(\omega_{k}-\omega)={\displaystyle\frac{\alpha\hbar}{2}}\omega_{c}^{1-\nu}\omega^{\nu}\Theta(\omega_{c}-\omega)$,
which is proportional to the dimensionless spin-bath coupling constant
$\alpha$. Unless otherwise stated, we take $\hbar=1$, $\nu=1$ (Ohmic case),
$\Delta=1$, cutoff energy $\omega_{c}=10$.

The QMC approach. The QMC method consists of quantum-to-classical mapping with
the introduction of an additional imaginary-time dimension, exact integration
of boson bath degrees of freedom and MC simulation of the resulting system
spin action up to the thermodynamic limit [1].

For the first step, we use the Suzuki-Trotter approximation [3] writing the
partition function as $Tr\left(e^{-\beta
H}\right)=\sum_{\\{\phi_{1},\phi_{N}\\}}\langle\phi_{1}|e^{-\frac{\beta}{N}H}|\phi_{2}\rangle\cdots\langle\phi_{N}|e^{-\frac{\beta}{N}H}|\phi_{1}\rangle,$
where $|\phi_{j}\rangle$ is the state of the system (both spins and bosons
degrees of freedom) at the $j$-th imaginary time. We therefore obtain a
classical system in $(1+1)$ dimensions, where $S_{i,j}=\pm 1$ is the value of
the spin at site $i=1\ldots L$ and time $j=1\ldots N$. For each site $i$ and
pair of imaginary times $j,j^{\prime}$, we introduce an auxiliary variable
$b_{i,j,j^{\prime}}=0,1$ that is equal to 1 if a phonon is emitted and
absorbed at $j$ and $j^{\prime}$, 0 otherwise.

After summing over the phonon degrees of freedom, the weight of a
configuration $\\{S_{i,j},b_{i,j,j^{\prime}}\\}$ is given by

$\displaystyle W(\\{S_{i,j},b_{i,j,j^{\prime}}\\})$
$\displaystyle=\exp\left[-{\cal H}_{nn}(\\{S_{i,j}\\})\right]\times$
$\displaystyle\prod_{i,j}\left(\delta_{B_{i,j},0}+\delta_{S_{i,j},-S_{i,j+1}}\delta_{B_{i,j},1}\right)\times$
$\displaystyle\prod_{j^{\prime}>j}\left[\frac{4}{\Delta^{2}}K_{S_{i,j},S_{i,j^{\prime}}}\left(\frac{\beta}{N}|j^{\prime}-j|\right)\right]^{b_{i,j,j^{\prime}}},$

where ${\cal H}_{nn}(\\{S_{i,j}\\})$ is a nearest neighbor classical
Hamiltonian induced by the quantum Hamiltonian (1),

${\cal
H}_{nn}(\\{S_{i,j}\\})=\sum_{i=1}^{L}\sum_{j=1}^{N}\left(-J_{\tau}S_{i,j}S_{i,j+1}-\tilde{J}S_{i,j}S_{i+1,j}\right),$
(6)

with $J_{\tau}=-\frac{1}{2}\ln\left(\frac{\tau\Delta}{2}\right)$,
$\tilde{J}=\frac{\tau J}{4}$, couplings along time and spatial direction,
respectively, and periodic boundary conditions in the time direction. The
variables $B_{i,j}$ are defined as
$B_{i,j}=\sum_{j^{\prime}}b_{i,j,j^{\prime}}$, and are limited to values 0 and
1 in the limit $N\to\infty$, since larger values are suppressed by factors
$\beta/N$. The long range kernel $K_{s,s^{\prime}}(\tau)$ is given by

$K_{s,s^{\prime}}(\tau)=\left\\{\begin{array}[]{ll}\gamma^{2}{\tilde{K}}_{s}(\tau)+(1-\gamma)^{2}{\tilde{K}}_{s^{\prime}}(\tau),&\mbox{if
$s\neq s^{\prime}$,}\\\
\rule{0.0pt}{15.0694pt}\gamma(1-\gamma)\left[{\tilde{K}}_{1}(\tau)+{\tilde{K}}_{-1}(\tau)\right],&\mbox{if
$s=s^{\prime}$,}\end{array}\right.$

with

${\tilde{K}}_{s}(\tau)=\int\limits_{0}^{\infty}\\!d\omega\,F(\omega)\frac{e^{s\omega\left(\tau-\frac{\beta}{2}\right)}}{e^{\omega\beta/2}-e^{-\omega\beta/2}}.$

Note that, for $\gamma\neq\frac{1}{2}$, the kernel $K_{s,s^{\prime}}$ breaks
the $Z_{2}$ symmetry, since $K_{1,-1}(\tau)\neq K_{-1,1}(\tau)$. In
particular, for $\gamma>\frac{1}{2}$, configurations with negative
magnetization $m_{z}$ are favoured, while for $\gamma<\frac{1}{2}$ it is the
opposite. Finally, in the classical limit, the function
${\tilde{K}}_{s}(\tau)$ becomes constant and independent of $s$:
${\tilde{K}}_{s}(\tau)=\frac{\alpha\omega_{c}}{2\beta}$, so that the $Z_{2}$
symmetry is restored [26].

Figure 1: Left Panel: Phase diagram with J vs. $\alpha$ for $\gamma=0.5$ in
the case of both quantum and classical baths. FERRO stands for ferromagnetic,
PARA for paramagnetic phase. Right Panel: Distribution function $P$ as a
function of the magnetization $m_{Z}$ for different values of $J$ in the case
of quantum baths at $\gamma=0.5$ and $\alpha=0.3$ (critical $J_{c}=2.47$). QMC
simulation parameters for right panel: $L=10$, and $\beta=10$.

As shown in the left panel of Fig. 1, for $\alpha=0$, QMC results recover the
well-known critical value $J_{c}=2$ [3]. We start from $\gamma=\frac{1}{2}$.
The local spin operator in the spin-bath Hamiltonian (4) becomes
$\gamma\sigma_{i}^{-}+(1-\gamma)\sigma_{i}^{+}=\sigma_{i}^{x}/2$, therefore it
is along the direction $x$ as the transverse field $\Delta$ in Eq. (2).
Clearly, at $J=0$, only the paramagnetic phase is stable. For finite $J$, we
always find a second order transition from a paramagnetic to a ferromagnetic
phase as a function of $\alpha$. Indeed, the critical $J_{c}$ gets enhanced
with increasing $\alpha$ along a transition line quite sensitive to system
parameters. Moreover, as detailed in the Supplemental Material, the dynamical
critical exponent $z$ is always equal to $1$ with changing $\alpha$, in
analogy with QMC results in Ising chains with Ohmic bond dissipation [27].
Therefore, our results differ from those of QMC literature [16, 17, 18] where
$\sigma^{z}_{i}$ is the local spin operator mediating the interaction with
quantum local baths [25]. QMC results can be interpreted within the VMF
approach, introduced below, in terms of an enhancement, proportional to
$\alpha$, of an effective transverse magnetic field.

The left panel of Fig. 1 shows that, as expected, the classical baths are less
effective in favouring ferromagnetic correlations. For example, at
$\alpha=0.3$, $J_{c}=2.47$ for quantum baths, while $J_{c}=3.67$ for classical
baths. In particular, for the quantum baths case, in the right panel of Fig.
1, we plot the distribution function $P$ as a function of the magnetization
$m_{Z}$ along z-axis for values of J smaller and slightly larger than $J_{C}$.
Actually, as expected for systems with $Z_{2}$ symmetry, the function $P$ is
symmetric around zero and shows a change from a mono-modal to a bimodal
character with crossing $J_{C}$.

Figure 2: Distribution function $P$ as a function of the magnetization $m_{Z}$
for different values of $J$ in the case of quantum (left panel) and classical
baths (right panel) for $\alpha=0.3$ and $\gamma=0.51$. QMC simulation
parameters: $L=10$, and $\beta=100$.

More interesting QMC results are obtained for $\gamma\neq\frac{1}{2}$ where
quantum bath fluctuations break $Z_{2}$ symmetry. As shown in the left panel
of Fig. 2, for $\gamma=0.51$, $P$ is centred around negative values of the
magnetization $m_{Z}$. With increasing $J$, the distribution $P$ shifts
towards more negative values. Therefore, as soon as $\gamma$ is different from
$\frac{1}{2}$, the effects of the bath induce a finite magnetization along
z-axis. QMC results for $\gamma\neq\frac{1}{2}$ can be better understood
through the VMF approach, exposed below, in terms of a dissipation induced
magnetic field along the z-axis. On the other hand, interactions between spins
and classical baths preserve $Z_{2}$ symmetry, therefore, as shown in the
right panel of Fig. 2, the function $P$ bears some resemblance with the
distribution function shown in the right panel of Fig. 1 at
$\gamma=\frac{1}{2}$.

The values of the bath parameter $\gamma$ smaller than $\frac{1}{2}$ favour
positive values of the magnetization along z-axis. In fact, configurations
with opposite values of $m_{z}$ are obtained for values of $\gamma$ symmetric
around $\frac{1}{2}$. For example, at $\gamma=0.49$ and fixed $\alpha$, the
distribution function P can be obtained from that shown in the left panel of
Fig. 2 making only the transformation $m_{Z}\rightarrow-m_{Z}$. In order to
interpret in a more effective way the numerically exact QMC results, in the
following we will present a zero temperature VMF approach within the framework
of the spin polaron [13]. Finally, by using both QMC and VMF approaches, we
will highlight the discontinuous behaviour of the magnetization $m_{Z}$ as a
function of $\gamma$ with crossing the value $\frac{1}{2}$ in the case of many
spins.

VMF: spin polaron framework. The trial wavefunction is chosen to be the ground
state of a test Hamiltonian obtained by replacing the term describing the
coupling between the spins in the original Hamiltonian with an effective
magnetic field along $z$ axis. In other words the trial wavefunction is:
$\left|\psi\right\rangle=\prod_{i=1}^{L}\left|\psi_{i}\right\rangle$, where
$\left|\psi_{i}\right\rangle$ represents the wavefunction of single spin
interacting with its local bath. The problem is then traced back to the
solution of one spin polaron that, in turn, will be variationally addressed
for both classical and quantum baths.

For the single spin, we will use two approaches: the first one is based on the
adiabatic approximation, rigorously valid in the classical limit; the second
one is a variational exact diagonalization method able to include quantum, non
adiabatic contributions. We simplify the notation for $L=1$ in Eq. (1):
$\sigma^{x}_{i}=\sigma_{x}$, $\sigma^{z}_{i}=\sigma_{z}$,
$\sigma^{+}_{i}=\sigma_{+}$, $\sigma^{-}_{i}=\sigma_{-}$, and, consequently,
$a_{i,k}^{+}=a_{k}^{+}$, $a_{i,k}=a_{k}$.

In the classical limit the lattice polarization cannot follow the spin
oscillations and the wave function of the system can be factorized into a
product of normalized variational functions $\left|\varphi\right\rangle$ and
$\left|g\right\rangle$, depending on the spin and bosonic coordinates,
respectively:
$\left|\psi\right\rangle=\left|\varphi\right\rangle\left|g\right\rangle$. The
expectation value of the Hamiltonian (1) on the state
$\left|\varphi\right\rangle$ provides an effective Hamiltonian for the bath
whose ground wave function is a coherent state:

$\left|g\right\rangle=e^{\sum_{k}\left(\frac{f\lambda_{k}}{\omega_{q}}a_{k}-h.c.\right)}\left|0\right\rangle.$
(7)

Here $\left|0\right\rangle$ is the bosonic vacuum state and
$f=\left\langle\varphi\left|\left[(1-\gamma)\sigma_{+}+\gamma\sigma_{-}\right]\right|\varphi\right\rangle$.
At this stage, the bosonic state $\left|g\right\rangle$ can be used to obtain
an effective Hamiltonian, $H_{eff}$, for the spin. It is straightforward to
show that $H_{eff}=-\frac{\tilde{\Delta}}{2}\sigma_{x}+C$, where
$\tilde{\Delta}=1+\sum_{k}\frac{2f\lambda^{2}_{k}}{\omega_{k}}$ and
$C=\sum_{k}\frac{f^{2}\lambda^{2}_{k}}{\omega_{k}}$, i.e. the adiabatic
approximation leads only to an enhancement of the transverse magnetic field
along $x$ axis, which is, however, enough to interpret the QMC results at
$\gamma=\frac{1}{2}$, shown in Fig. 1.

Figure 3: Magnetization $m_{Z}$ as a function of $\gamma$ for different values
of $J$ and $\alpha$ by using both QMC and VMF methods. QMC simulation
parameters: $\beta=100$; $L=1$ for (a) and (b), $L=40$ for (c), and $L=10$ for
(d).

In order to get a proper inclusion of the non-adiabatic contributions,
relevant for $\gamma\neq\frac{1}{2}$, we first apply an unitary transformation
$H_{1}=e^{S_{1}}He^{-S_{1}}$, with
$S_{1}=-\sum_{k}\left(\frac{f\lambda_{k}}{\omega_{k}}a_{k}-h.c.\right)$, as
suggested by Eq. (7). The transformed Hamiltonian assumes the form:

$H_{1}=H_{eff}+\sum_{q}\omega_{k}a^{\dagger}_{k}a_{k}+\tilde{H}_{I},$ (8)

where $\tilde{H}_{I}=\sum_{k}\lambda_{k}Aa_{k}+h.c$, $A$ being an operator
acting only on the spin subsystem:
$A=\left[(1-\gamma)\sigma_{+}+\gamma\sigma_{-}\right]-f$. By treating
$\tilde{H}_{I}$ as perturbation, one includes the non adiabatic contributions.
The first order of the perturbation theory, followed by the assumption of no
correlation between the emission of different bosons, suggests the following
trial state for the bath: $e^{-S_{2}}\left|0\right\rangle$, where
$S_{2}=-\sum_{k}\sigma_{z}\left(\frac{f_{1}\lambda_{k}}{\omega_{k}+\tilde{\Delta}}a_{k}-h.c.\right)$
and $f_{1}=\gamma-1/2$. In other words, $S_{2}$ takes into account the
possibility that the bath can follow instantaneously the spin oscillations.
The effective Hamiltonian for the spin turns out to be:

$H_{eff}=-\frac{\Delta_{eff}}{2}\sigma_{x}+C_{eff}+h_{eff}\sigma_{z},$ (9)

where
$\Delta_{eff}=\left(\tilde{\Delta}+\sum_{k}\frac{4f_{1}^{2}\lambda^{2}_{k}}{\omega_{k}+\tilde{\Delta}}\right)e^{-\sum_{k}\frac{2f_{1}^{2}\lambda_{k}^{2}}{(\omega_{k}+\tilde{\Delta})^{2}}}$,
$C_{eff}=C+\sum_{k}\frac{\omega_{k}f_{1}^{2}\lambda^{2}_{k}}{\left(\omega_{k}+\tilde{\Delta}\right)^{2}}$,
and
$h_{eff}=2f\sum_{k}\frac{f_{1}\lambda^{2}_{k}}{\omega_{k}+\tilde{\Delta}}$. It
is evident that the quantum fluctuations of the bath induce an effective
magnetic field $h_{eff}$ along $z$ axis that breaks the $Z_{2}$ symmetry.
Moreover, the field $h_{eff}$ changes sign with crossing $\gamma=\frac{1}{2}$.

The VMF approach can be further improved by by diagonalizing the original
Hamiltonian in the subspace spanned by the following $2M$ wavefunctions:
$\psi_{i}=e^{\sum_{k}\left(l_{k,i}a_{k}-h.c.\right)}\left|0\right\rangle\left|\uparrow\right\rangle$
and
$\phi_{i}=e^{\sum_{k}\left(h_{k,i}a_{k}-h.c.\right)}\left|0\right\rangle\left|\downarrow\right\rangle$,
with $i=1,....M$, in analogy with approaches for the Holstein model [21, 22,
23, 24]. Here, as suggested by the above described approach, we assume:
$l_{k,i}=\frac{f\lambda_{k}}{\omega_{k}}+\frac{\lambda_{k}f_{i}}{\omega_{k}+\Delta_{i}}$,
$h_{k,i}=\frac{f\lambda_{k}}{\omega_{k}}-\frac{\lambda_{k}f_{i}}{\omega_{q}+\Delta_{i}}$,
being $f$, $f_{i}$ and $\Delta_{i}$ $2M+1$ variational parameters. In Fig.
3(a) and 3(b), where we successfully compare the calculated magnetization
$m_{Z}$ with QMC results for the single spin, we prove that $M=2$ is enough to
get a very accurate description of the spin polaron in any regime. Moreover,
in the case of a single spin, as expected, the magnetization vanishes at
$\gamma=\frac{1}{2}$, and, as a function of $\gamma$, shows a crossover
behaviour, weakly dependent on $\alpha$, from positive to negative values.

In the case of a spin chain, as shown in Fig. 3(c) and (d), the calculated
magnetization $m_{Z}$ perfectly matches QMC results for $J$ larger than
$J_{C}$ (panel (c), $\alpha=0.3$, $J_{C}=2.47$ at $\gamma=\frac{1}{2}$) and
smaller than $J_{C}$ (panel (d), $\alpha=1$, $J_{C}=3.63$ at
$\gamma=\frac{1}{2}$). In the first case, there is an actual first order
transition induced by $\gamma$ between states with opposite magnetization,
while, in the second case, only a crossover takes place at
$\gamma=\frac{1}{2}$. We remark that the first order transition occurs also
for lower values of $\alpha$, therefore the spin system shows a large
sensitivity to external perturbations even at minimal coupling [28].

Figure 4: Phase diagram with $\gamma$ vs. $\alpha$ at $J=2.5$ (left panel),
with $\gamma$ vs. J at $\alpha=0.3$ (right panel). Vertical arrows indicate
first order transitions, circles denote second order critical points which can
be encircled (curved arrows). QMC simulation parameters: $\beta=100$, $L=40$.

Finally, we complement the phase diagram of Fig. 1, showing in Fig. 4 the
diagram with $\gamma$ vs. $\alpha$ at fixed $J$ (left panel ), and the diagram
with $\gamma$ vs. $J$ at fixed $\alpha$ (right panel). It is apparent that one
gets an horizontal line at $\gamma=\frac{1}{2}$ of first order transitions
terminating with a second order critical point (denoted with a circle in the
panels of Fig. 4). From the comparison with the phase diagram in Fig. 1, the
critical point at $\alpha=0.3$ in the right panel corresponds to $J_{c}=2.47$.
As expected in this kind of transition [1], one can go continuously from a
magnetization state to its opposite encircling the second order critical point
(curved arrows).

Conclusions. In this letter, through a numerically exact QMC method and VMF
approach, we have shown that a combination of Jaynes-Cummings and anti-Jaynes-
Cummings spin-bath couplings drive the transition in a transverse field Ising
chain from the second to the first order even for a low dissipation strength.
These features are present only for quantum baths. We remark that, in addition
to spin-spin interactions [9, 10], Jaynes-Cummings couplings can be
implemented in quantum simulators [29, 30], which, therefore, can potentially
test the results exposed in this work. Moreover, the spin system is not only
sensitive to the properties of non-classical baths, but also to weak magnetic
fields, implying that the set-up proposed in this work can be also relevant
for the realization of small-scale quantum sensors. Finally, the approaches
used in this work can be easily generalized to more complex spin systems [31,
32] and spin-bath couplings [33], and in principle allow also to include the
effects of coherent sinusoidal drives [34].

###### Acknowledgements.

V.C. and G.D.F. acknowledge a critical reading of the manuscript by N.
Nagaosa.

## References

  * [1] S. Sachdev, Quantum Phase Transitions 2nd edn, Cambridge University Press (2011).
  * [2] D. Rossini, and E. Vicari, Phys. Rep. 96, 1 (2021).
  * [3] S. Suzuki, J. Inoue, and B. K. Chakrabarti, Quantum Ising Phases and Transitions in Transverse Ising Models 2nd edn, Springer (2012).
  * [4] L. M. Sieberer, M. Buchhold, and S. Diehl, Rep. Prog. Phys. 79, 096001 (2016).
  * [5] H.-P. Breuer, and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, New York (2007) .
  * [6] A. O. Caldeira, and A. J. Leggett, Ann. Phys. 149, 374 (1983).
  * [7] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * [8] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * [9] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletic, and M. D. Lukin, Nature (London) 551, 579 (2017).
  * [10] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletic, and M. D. Lukin, Nature 568, 207 (2019).
  * [11] S. Ebadi1, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletic, and M. D. Lukin, Nature 595, 227 (2021).
  * [12] C. Noh and D. G. Angelakis, Rep. Prog. Phys. 80, 016401 (2017).
  * [13] K. Le Hur, in L. D. Carr (Ed.), Understanding Quantum Phase Transitions, CRC Press (2010).
  * [14] G. De Filippis, A. de Candia, A. S. Mishchenko, L. M. Cangemi, A. Nocera, P. A. Mishchenko, M. Sassetti, R. Fazio, N. Nagaosa, and V. Cataudella, Phys. Rev. B 104, 060410 (2021).
  * [15] G. De Filippis, A. de Candia, G. Di Bello, C. A. Perroni, L. M. Cangemi, A. Nocera, M. Sassetti, R. Fazio, and V. Cataudella, arXiv:2205.11555 (2022).
  * [16] P. Werner, M. Troyer, and S. Sachdev, J. Phys. Soc. Jpn. Suppl. 74, 67 (2005).
  * [17] P. Werner, K. Volker, M. Troyer, and S. Chakravarty, Phys. Rev. Lett. 94, 047201 (2005).
  * [18] P. Werner, and M. Troyer, Prog. Theor. Phys. Suppl. 160, 395 (2005).
  * [19] G. M. Uhdre, D. Cius, and F. M. Andrade, Phys. Rev. A 105, 013703 (2022).
  * [20] J. A. Omolo, arXiv:2103.09546 (2021).
  * [21] G. De Filippis, V. Cataudella, V. Marigliano Ramaglia, and C. A. Perroni, Phys. Rev. B 72, 014307 (2005).
  * [22] G. De Filippis, V. Cataudella, A. S. Mishchenko, C. A. Perroni, and N. Nagaosa, Phys. Rev. B 80, 195104 (2009).
  * [23] A. Nocera, M. Soltanieh-ha, C. A. Perroni, V. Cataudella, and A. E. Feiguin, Phys. Rev. B 90, 195134 (2014).
  * [24] C. A. Perroni, D. Ninno, and V. Cataudella, J. Phys. Condens. Matter 28, 373001 (2016).
  * [25] We have checked that our method perfectly reproduces the QMC results reported in [17].
  * [26] In the Supplemental Material, we will provide additional details about the QMC method, expanding not only the limit of classical bath oscillators, but introducing also Binder cumulants, which are used to evaluate the phase diagram shown in Fig. 1.
  * [27] I. B. Sperstad, E. B. Stiansen, and A. Sudbo, Phys. Rev. B 81, 104302 (2010).
  * [28] In the Supplemental Material, we will discuss the finite size scaling of the magnetization $m_{Z}$ in the parameter regime of Fig. 3(c) and the interesting behaviour of the magnetization $m_{X}$ for the same parameters as in Fig. 3.
  * [29] J. S. Pedernales, I. Lizuain, S. Felicetti, G. Romero, L. Lamata, and E. Solano, Sci. Rep. 5, 15472 (2015).
  * [30] J. Braumuller, M. Marthaler, A. Schneider, A. Stehli, H. Rotzinger, M. Weides, and A. V. Ustinov, Nat. Commun. 8, 779 (2017).
  * [31] J. Jin, W.-B. He, F. Iemini, D. Ferreira, Y.-D. Wang, S. Chesi, and R. Fazio, Phys. Rev. B 104, 214301 (2021).
  * [32] M. Weber, D. J. Luitz, and F. F. Assaad, arXiv:2112.02124v1 (2021).
  * [33] Y. Zheng and S. Yang, Phys. Rev. A 104, 023304 (2021).
  * [34] H. Landa, M. Schiró, and G. Misguich, Phys. Rev. Lett. 124, 043601 (2020).