The intensity $\tilde{I}_{\nu,i}$ emerging from the surface is the resulting
radiation as having propagated from the interior. This radiation originates
from all points beneath the surface as attenuated on a scale set by the
electric skin depth $\delta_{{\rm elec}}$, where

$\delta_{{\rm elec}}\equiv\frac{\lambda}{4\pi\kappa},$ (D4)

in which $\kappa={\rm Im}(n)$ and where $n$ is the complex index of refraction
of the medium that is related to the dielectric constant $\epsilon$ via
$\epsilon\equiv n^{2}$, and is generally a function of a photon’s wavelength.
While the electric skin depth is temperature dependent, for simplicity
hereafter we will treat it as insensitive to temperature. A material’s
$\epsilon$ is usually measured in the laboratory (e.g., for cometary ices at
radar wavelengths see Heggy et al., 2012) and is tabulated in terms of its
real and imaginary parts for which we use the convention in de Kleer et al.
(2021), $\varepsilon^{\prime}\equiv{\rm Re}(\epsilon)$ and
$\varepsilon^{\prime\prime}\equiv{\rm Im}(\epsilon)$. In all of our
calculations below, we assume that $\epsilon$ has no subsurface variation. As
such

$\tilde{I}_{\nu,i}(\cos\theta_{i})=\frac{\displaystyle\left<E_{{\rm
eff}}\right>\int_{0}^{\infty}I_{\nu,i}(z)e^{-z\big{/}\left(\delta_{{\rm
elec}}\cos\theta_{t,i}\right)}dz}{\displaystyle\int_{0}^{\infty}e^{-z\big{/}\left(\delta_{{\rm
elec}}\cos\theta_{t,i}\right)}dz}$ (D5)

where $z$ is depth as measured normal to facet $i$’s surface, $\theta_{t,i}$
is the angle formed between a subsurface propagating photon’s path and facet
$i$’s unit normal (for reference see Fig. 23). We have formally included an
X-band emissivity factor $\left<E_{{\rm eff}}\right>$ as is common practice
when dealing with solar system ices of unknown structural and transmission
properties. Because we have assumed $\epsilon$ is constant in the subsurface
the photon’s subsurface slant-path trajectory changes upon exiting into vacuum
and, therefore, $\theta_{t,i}$ relates to $\theta_{i}$ via Snell’s Law, which
following its reformulation in de Kleer et al. (2021), is given by

$\cos\theta_{t,i}=\sqrt{1-\frac{\sin^{2}\theta_{i}}{\varepsilon^{\prime}}}.$
(D6)

The subsurface intensity $I_{\nu,i}$ is

$I_{\nu,i}(z)=\frac{\displaystyle
2h\nu^{3}\big{/}c^{2}}{\displaystyle{\exp\left[{{\frac{h\nu}{kT_{i}(z)}}}\right]-1}}\approx\frac{2k\nu^{2}}{c^{2}}T_{i}(z),$
(D7)

The Rayleigh-Jeans limit is the approximate form given on the RHS of Eq. (D7).
The subsurface temperature profile beneath facet $i$, $T_{i}=T_{i}(z,{\cal
I},\cdots)$, are those as determined in previous sections, and we note that
these temperature solutions depends upon the material’s thermal inertia
(${\cal I}$), which we shall keep general for our considerations forthwith.

Using the Rayleigh-Jeans limiting form of Eq. (D7) in (D5) permits the re-
expression of Eq. (D3) in terms of a model prediction brightness temperature,
hereafter simply as $T_{b}$,

$\displaystyle\tilde{I}_{\nu}$ $\displaystyle=$
$\displaystyle\frac{2k\nu^{2}}{c^{2}}T_{b}\big{(}{\cal
I},\varepsilon^{\prime},\varepsilon^{\prime\prime},\big{<}E_{\rm
eff}\big{>}\big{)},$ (D8)

where

$T_{b}=T_{b}\Big{(}{\cal
I},\varepsilon^{\prime},\varepsilon^{\prime\prime},\big{<}E_{\rm
eff}\big{>}\Big{)}=\overline{T}_{b,i}\equiv\frac{\displaystyle{\sum_{\forall
i\ {\rm visible}}\sigma_{i}T_{b,i}}}{\displaystyle{\sum_{\forall i\ {\rm
visible}}\sigma_{i}}},$ (D9)

with a predicted X-band brightness temperature for each visible facet $i$,
$T_{b,i}$ defined by

$T_{b,i}\Big{(}{\cal
I},\varepsilon^{\prime},\varepsilon^{\prime\prime},\big{<}E_{\rm
eff}\big{>}\Big{)}\equiv\frac{\displaystyle{\big{<}E_{\rm
eff}\big{>}\int_{0}^{\infty}T_{i}(z,{\cal I})e^{-z\big{/}\left(\delta_{{\rm
elec}}\cos\theta_{t,i}\right)}}dz}{\delta_{{\rm elec}}\cos\theta_{t,i}}.$
(D10)

The overline symbol is intended here as a short-hand representation of the
facet-area weighted averaging. For any given instantaneous temperature
solution determined using the methods of section 3, we can produce a model
prediction $T_{b}$ that can be compared directly against the REX observed
value $T_{b,{\rm obs}}$. $T_{b}$ is a function of four unknown parameters,
i.e., ${\cal I}$ (through $T_{i}$), the real and imaginary parts of the
dielectric constant $\epsilon$, and the X-band emissivity $\big{<}E_{\rm
eff}\big{>}$. The dielectric constant determines the electric skin depth
($\delta_{{\rm elec}}$) and the photon’s subsurface propagation angle
$\theta_{t,i}$ that leads to a photon emission angle $\theta_{i}$ that, upon
leaving Arrokoth’s surface, leads directly to the spacecraft. Each
$\theta_{t,i}$ is calculated by Snell’s Law, Eq. (D6), which is a known
function of the set of viewing angles $\theta_{i}$ formed between each facet’s
unit normal ${\bf{\hat{n}}}_{i}$ and the line connecting the spacecraft’s
position to each facet’s center.

We note that using an averaging approach Bird et al. (2022) develop $T_{b}$
along a parallel route that short-circuits the large scale computations
involved in developing Eq. (D9) from Eq. (D10). Instead, $T_{b}$ is estimated
straight away from (D10) where $T_{i}$ is replaced by a visible facet surface
area weighted temperature, $\overline{T(z)}$, and $\theta_{t,i}$ is replaced
with a similarly facet area weighted average of the spacecraft-surface-normal
inclination angle, $\overline{\theta_{t}}$. In other words

$T_{b}\approx\frac{\displaystyle{\big{<}E_{\rm
eff}\big{>}\int_{0}^{\infty}\overline{T_{i}(z,{\cal
I})}e^{-z\big{/}\left(\delta_{{\rm
elec}}\overline{\cos\theta_{t,i}}\right)}}dz}{\delta_{{\rm
elec}}\overline{\cos\theta_{t,i}}}.$ (D11)

Formally speaking, the two approaches are not-commensurate as the average of
products – namely that of the temperature and the depth dependent exponential
factor – is not necessarily equal to the product of their averages. However,
we do find that in this application the final results are largely consistent
with one another.

Thus, Following the procedures described in Section 3, we construct a set of
model temperature solutions for several values of ${\cal I}_{j}$ contained in
the range of $0.5$tiu and $300$tiu. We collect the set of predicted
temperature solutions for encounter day $T_{i}(z,{\cal I}_{j})$ and generate
predicted values of $T_{b}$ for the range of dielectric values
$1<\varepsilon^{\prime}<3.5$ and $10^{-4}<\epsilon^{\prime\prime}<10$, which
are motivated by those values of $\epsilon$ thought to be relevant to cometery
ices. The numerical solutions determined for each $T_{i}(z,{\cal I}_{j})$ are
tabulated only down to a depth $z_{b}=10$m because we find that the
temperatures below that depth have asymptoted to their orbitally averaged
values, i.e., $T_{i}(z>z_{b},{\cal I}_{j})\approx T_{i}(z_{b},{\cal I}_{j})$.
As such, the integral in Eq. (D10) is well-approximated with

$\displaystyle\frac{\displaystyle{\int_{0}^{\infty}T_{i}(z,{\cal
I}_{j})e^{-z\big{/}\left(\delta_{{\rm
elec}}\cos\theta_{t,i}\right)}dz}}{\delta_{{\rm
elec}}\cos\theta_{t,i}}\approx\frac{\displaystyle{\int_{0}^{z_{b}}T_{i}(z,{\cal
I}_{j})e^{-z\big{/}\left(\delta_{{\rm
elec}}\cos\theta_{t,i}\right)}dz}}{\delta_{{\rm
elec}}\cos\theta_{t,i}}+T_{i}(z_{b},{\cal
I}_{j})e^{-z_{b}\big{/}\left(\delta_{{\rm elec}}\cos\theta_{t,i}\right)}.$
(D12)

We observe that the added correction term on the RHS of the above expression
is important only when $\delta_{\rm{elec}}\cos\theta_{t,i}\ll z_{b}$.

Fig. 15 displays predicted $T_{b}$ for eight selected values of ${\cal I}$.
Each color contour plot shows $T_{b}$ as a function of $\varepsilon^{\prime}$
and $\varepsilon^{\prime\prime}$. Drawn atop each are black contour lines
representing $T_{b,{\rm obs}}=29$K, with dashed black contours showing the
error bounds at 34K and 24K. Also drawn with white contours a level sets for
four corresponding values of $\delta_{{\rm elec}}=10^{-2},10^{-1},1,10$m,
nominally indicating corresponding depths from which the bulk of the radiation
emerges for a given value of $\varepsilon$ of interest. The qualities and
implications of these solutions are discussed further in the body of the text.
However we do observe from Fig. 15 that the predicted values of $T_{b}$ are
weakly dependent on $\varepsilon^{\prime}$ in the range 1 and 3.5.

## References

  * Andersson & Suga (1994) Andersson, O., & Suga, H. 1994, Solid State Communications, 91, 985, doi: 10.1016/0038-1098(94)90438-3
  * Attree et al. (2018) Attree, N., Groussin, O., Jorda, L., et al. 2018, A&A, 610, A76, doi: 10.1051/0004-6361/201731937
  * Beyer et al. (2019a) Beyer, R., Porter, S., Schenk, P., et al. 2019a, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–849
  * Beyer et al. (2019b) Beyer, R. A., Weaver, H. A., Porter, S. B., et al. 2019b, in Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 2258
  * Bierson & Nimmo (2019) Bierson, C. J., & Nimmo, F. 2019, Icarus, 326, 10, doi: 10.1016/j.icarus.2019.01.027
  * Brouet et al. (2015) Brouet, Y., Levasseur-Regourd, A. C., Sabouroux, P., et al. 2015, A&A, 583, A39, doi: 10.1051/0004-6361/201526099
  * Buie et al. (2020) Buie, M. W., Porter, S. B., Tamblyn, P., et al. 2020, AJ, 159, 130, doi: 10.3847/1538-3881/ab6ced
  * Castillo-Rogez et al. (2019) Castillo-Rogez, J., Vernazza, P., & Walsh, K. 2019, MNRAS, 486, 538, doi: 10.1093/mnras/stz786
  * Castillo-Rogez et al. (2012) Castillo-Rogez, J. C., Johnson, T. V., Thomas, P. C., et al. 2012, Icarus, 219, 86, doi: 10.1016/j.icarus.2012.02.002
  * Cruikshank et al. (2005a) Cruikshank, D. P., Imanaka, H., & Dalle Ore, C. M. 2005a, Advances in Space Research, 36, 178, doi: 10.1016/j.asr.2005.07.026
  * Cruikshank et al. (2005b) Cruikshank, D. P., Owen, T. C., Dalle Ore, C., et al. 2005b, Icarus, 175, 268, doi: 10.1016/j.icarus.2004.09.003
  * de Kleer et al. (2021) de Kleer, K., Butler, B., de Pater, I., et al. 2021, , 2, 5, doi: 10.3847/PSJ/abcbf4
  * Earle et al. (2017) Earle, A. M., Binzel, R. P., Young, L. A., et al. 2017, Icarus, 287, 37, doi: 10.1016/j.icarus.2016.09.036
  * Efimov et al. (2011) Efimov, V. B., Izotov, A. N., Levchenko, A. A., Mezhov-Deglin, L. P., & Khasanov, S. S. 2011, JETP Letters, 94, 621
  * Estrada et al. (2016) Estrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200, doi: 10.3847/0004-637X/818/2/200
  * Ferrari (2018) Ferrari, C. 2018, Space Sci. Rev., 214, 111, doi: 10.1007/s11214-018-0546-x
  * Ferrari & Lucas (2016) Ferrari, C., & Lucas, A. 2016, A&A, 588, A133, doi: 10.1051/0004-6361/201527625
  * Ferrari et al. (2021) Ferrari, C., Lucas, A., & Jacquemoud, S. 2021, A&A, 655, A8, doi: 10.1051/0004-6361/202141223
  * Flasar et al. (2005) Flasar, F. M., Achterberg, R. K., Conrath, B. J., et al. 2005, Science, 307, 1247, doi: 10.1126/science.1105806
  * Fray & Schmitt (2009) Fray, N., & Schmitt, B. 2009, Planet. Space Sci., 57, 2053, doi: 10.1016/j.pss.2009.09.011
  * Gläser & Gläser (2019) Gläser, P., & Gläser, D. 2019, A&A, 627, A129, doi: 10.1051/0004-6361/201935514
  * Grundy et al. (2016) Grundy, W. M., Binzel, R. P., Buratti, B. J., et al. 2016, Science, 351, aad9189, doi: 10.1126/science.aad9189
  * Grundy et al. (2020) Grundy, W. M., Bird, M. K., Britt, D. T., et al. 2020, Science, 367, aay3705, doi: 10.1126/science.aay3705
  * Gundlach & Blum (2012) Gundlach, B., & Blum, J. 2012, Icarus, 219, 618, doi: 10.1016/j.icarus.2012.03.013
  * Gusarov et al. (2003) Gusarov, A., Laoui, T., Froyen, L., & Titov, V. 2003, Int. J. Heat Mass Transfer, 46, 1103, doi: 10.1016/S0017-9310(02)00370-8
  * Hartmann (1987) Hartmann, W. K. 1987, Icarus, 71, 57, doi: 10.1016/0019-1035(87)90162-X
  * Heggy et al. (2012) Heggy, E., Palmer, E. M., Kofman, W., et al. 2012, Icarus, 221, 925, doi: 10.1016/j.icarus.2012.09.023
  * Hofgartner et al. (2021) Hofgartner, J. D., Buratti, B. J., Benecchi, S. D., et al. 2021, Icarus, 356, 113723, doi: 10.1016/j.icarus.2020.113723
  * Howett et al. (2010) Howett, C. J. A., Spencer, J. R., Pearl, J., & Segura, M. 2010, Icarus, 206, 573, doi: 10.1016/j.icarus.2009.07.016
  * Hu et al. (2017) Hu, X., Shi, X., Sierks, H., et al. 2017, MNRAS, 469, S295, doi: 10.1093/mnras/stx1607
  * Jackson & Desch (2021) Jackson, A. P., & Desch, S. J. 2021, Journal of Geophysical Research (Planets), 126, e06706, doi: 10.1029/2020JE006706
  * Jagota et al. (2012) Jagota, A., Paretkar, D., & Ghatak, A. 2012, Phys. Rev. E, 85, 051602, doi: 10.1103/PhysRevE.85.051602
  * Johnson et al. (1971) Johnson, K. L., Kendall, K., & Roberts, A. D. 1971, Proceedings of the Royal Society of London Series A, 324, 301, doi: 10.1098/rspa.1971.0141
  * Johnson & Lunine (2005) Johnson, T. V., & Lunine, J. I. 2005, Nature, 435, 69, doi: 10.1038/nature03384
  * Keane et al. (2020) Keane, J. T., Umurhan, O. M., Porter, S. B., et al. 2020, in Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 2444
  * Klinger (1980) Klinger, J. 1980, Science, 209, 271, doi: 10.1126/science.209.4453.271
  * Le Gall et al. (2014) Le Gall, A., Leyrat, C., Janssen, M. A., et al. 2014, Icarus, 241, 221, doi: 10.1016/j.icarus.2014.06.011
  * Le Gall et al. (2016) Le Gall, A., Malaska, M. J., Lorenz, R. D., et al. 2016, Journal of Geophysical Research (Planets), 121, 233, doi: 10.1002/2015JE004920
  * Lebofsky (1975) Lebofsky, L. A. 1975, Icarus, 25, 205, doi: 10.1016/0019-1035(75)90020-2
  * Lellouch et al. (2013) Lellouch, E., Santos-Sanz, P., Lacerda, P., et al. 2013, A&A, 557, A60, doi: 10.1051/0004-6361/201322047
  * Lisse et al. (2021) Lisse, C. M., Young, L. A., Cruikshank, D. P., et al. 2021, Icarus, 356, 114072, doi: 10.1016/j.icarus.2020.114072
  * Matson et al. (2009) Matson, D. L., Castillo-Rogez, J. C., Schubert, G., Sotin, C., & McKinnon, W. B. 2009, The Thermal Evolution and Internal Structure of Saturn’s Mid-Sized Icy Satellites, ed. M. K. Dougherty, L. W. Esposito, & S. M. Krimigis, 577, doi: 10.1007/978-1-4020-9217-6_18
  * McKinnon et al. (2020) McKinnon, W. B., Richardson, D. C., Marohnic, J. C., et al. 2020, Science, 367, aay6620, doi: 10.1126/science.aay6620
  * Mellon (1997) Mellon, M. T. 1997, J. Geophys. Res., 102, 25617, doi: 10.1029/97JE02582
  * Meyer et al. (2012) Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2012, MNRAS, 422, 2102, doi: 10.1111/j.1365-2966.2012.20744.x
  * Molaro et al. (2015) Molaro, J. L., Byrne, S., & Langer, S. A. 2015, Journal of Geophysical Research (Planets), 120, 255, doi: 10.1002/2014JE004729
  * Müller et al. (2020) Müller, T., Lellouch, E., & Fornasier, S. 2020, in The Trans-Neptunian Solar System, ed. D. Prialnik, M. A. Barucci, & L. Young, 153–181, doi: 10.1016/B978-0-12-816490-7.00007-2
  * Ostro et al. (2006) Ostro, S. J., West, R. D., Janssen, M. A., et al. 2006, Icarus, 183, 479, doi: 10.1016/j.icarus.2006.02.019
  * Paillou et al. (2008) Paillou, P., Lunine, J., Ruffié, G., et al. 2008, Geophys. Res. Lett., 35, L18202, doi: 10.1029/2008GL035216
  * Pätzold et al. (2016) Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63, doi: 10.1038/nature16535
  * Piqueux & Christensen (2009) Piqueux, S., & Christensen, P. R. 2009, Journal of Geophysical Research (Planets), 114, E09005, doi: 10.1029/2008JE003308
  * Porter et al. (2019) Porter, S., Beyer, R., Keane, J., et al. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–311
  * Porter et al. (2018) Porter, S. B., Buie, M. W., Parker, A. H., et al. 2018, AJ, 156, 20, doi: 10.3847/1538-3881/aac2e1
  * Richardson & Bowling (2014) Richardson, J. E., & Bowling, T. J. 2014, Icarus, 234, 53, doi: 10.1016/j.icarus.2014.02.015
  * Richardson et al. (2019) Richardson, J. E., Graves, K. J., Harris, A. W., & Bowling, T. J. 2019, Icarus, 329, 207, doi: 10.1016/j.icarus.2019.03.027
  * Rivera-Valentin et al. (2011) Rivera-Valentin, E. G., Blackburn, D. G., & Ulrich, R. 2011, Icarus, 216, 347, doi: 10.1016/j.icarus.2011.09.006
  * Rozitis & Green (2011) Rozitis, B., & Green, S. F. 2011, MNRAS, 415, 2042, doi: 10.1111/j.1365-2966.2011.18718.x
  * Schloerb et al. (2015) Schloerb, F. P., Keihm, S., von Allmen, P., et al. 2015, A&A, 583, A29, doi: 10.1051/0004-6361/201526152
  * Schmitt et al. (2017) Schmitt, B., Philippe, S., Grundy, W. M., et al. 2017, Icarus, 287, 229, doi: 10.1016/j.icarus.2016.12.025
  * Shoshany et al. (2002) Shoshany, Y., Prialnik, D., & Podolak, M. 2002, Icarus, 157, 219, doi: 10.1006/icar.2002.6815
  * Shulman (2004) Shulman, L. M. 2004, A&A, 416, 187, doi: 10.1051/0004-6361:20031746
  * Simon et al. (2018) Simon, J. I., Cuzzi, J. N., McCain, K. A., et al. 2018, Earth and Planetary Science Letters, 494, 69, doi: 10.1016/j.epsl.2018.04.021
  * Spencer et al. (2004) Spencer, J. R., Pearl, J. C., Segura, M., & Cassini CIRS. 2004, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 36, AAS/Division for Planetary Sciences Meeting Abstracts #36, 04.07
  * Spencer et al. (2020) Spencer, J. R., Stern, S. A., Moore, J. M., et al. 2020, Science, 367, aay3999, doi: 10.1126/science.aay3999
  * Stern et al. (2021) Stern, S. A., Keeney, B., Singer, K. N., et al. 2021, , 2, 87, doi: 10.3847/PSJ/abee26
  * Stern et al. (2019) Stern, S. A., Weaver, H. A., Spencer, J. R., et al. 2019, Science, 364, aaw9771, doi: 10.1126/science.aaw9771
  * Style et al. (2013) Style, R. W., Boltyanskiy, R., Che, Y., et al. 2013, Phys. Rev. Lett., 110, 066103, doi: 10.1103/PhysRevLett.110.066103
  * Titus & Cushing (2012) Titus, T. N., & Cushing, G. E. 2012, in Third International Planetary Dunes Workshop: Remote Sensing and Data Analysis of Planetary Dunes, Vol. 1673, 95–96
  * Tosi et al. (2019) Tosi, F., Capaccioni, F., Capria, M. T., et al. 2019, Nature Astronomy, 3, 649, doi: 10.1038/s41550-019-0740-0
  * Umurhan et al. (2019) Umurhan, O. M., Keane, J. T., Porter, S. B., et al. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–749
  * Verbiscer et al. (2009) Verbiscer, A. J., Skrutskie, M. F., & Hamilton, D. P. 2009, Nature, 461, 1098, doi: 10.1038/nature08515
  * White et al. (2016) White, O. L., Umurhan, O. M., Moore, J. M., & Howard, A. D. 2016, Journal of Geophysical Research (Planets), 121, 21, doi: 10.1002/2015JE004846


# The perils of being unhinged:
On the accuracy of classifiers minimizing
a noise-robust convex loss

Philip M. Long
Google
<EMAIL_ADDRESS>
Rocco A. Servedio
Columbia University
<EMAIL_ADDRESS>

###### Abstract

van Rooyen et al. (2015) introduced a notion of convex loss functions being
robust to random classification noise, and established that the “unhinged”
loss function is robust in this sense. In this note we study the accuracy of
binary classifiers obtained by minimizing the unhinged loss, and observe that
even for simple linearly separable data distributions, minimizing the unhinged
loss may only yield a binary classifier with accuracy no better than random
guessing.

##  1 Introduction

As van Rooyen et al. noted in the first sentence of the abstract of van Rooyen
et al. (2015), “Convex potential minimisation is the _de facto_ approach to
binary classification.” Given the ubiquity of this approach, it is natural to
study its abilities and limitations in the presence of noise, and indeed this
is the subject of many works (see Zhang, 2004; Bartlett et al., 2006; Long and
Servedio, 2010; Manwani and Sastry, 2013; Natarajan et al., 2013; van Rooyen
et al., 2015; Ghosh et al., 2017)).

The aim of this note is to clarify the connection between minimizing a convex
potential function which is “robust to classification noise” in the sense of
van Rooyen et al. (2015), and learning (i.e. performing accurate
classification).

Background. Motivated by the observation that the popular AdaBoost algorithm
(which works by minimizing the (convex) exponential potential function) can
have empirically poor classification accuracy when run on noisy data
(Dietterich, 2000; Freund and Schapire, 1996; Maclin and Opitz, 1997), Long
and Servedio (2010) studied the performance of classification algorithms which
work by minimizing a convex potential function in settings where linearly
separable data is contaminated with random classification noise (RCN). The
main result of Long and Servedio (2010) is a proof that for a certain simple
learning problem corresponding to a “clean” data distribution ${\cal D}_{1}$
that is linearly separable with a margin, for _any_ “convex potential
function” $\phi$, minimizing $\phi$ over all linear combinations of base
features in the presence of random classification noise only yields a binary
classifier with an error rate of 1/2 under the clean distribution ${\cal
D}_{1}$. (Here a “convex potential function” is a convex function
$\phi:\mathds{R}\to\mathds{R}$ satisfying certain mild conditions which we
detail in Definition 1.) This is in sharp contrast with the fact that, in the
noise-free setting of a data distribution that is linearly separable with a
margin, driving the potential to zero leads to a perfectly accurate binary
classifier.

In an effort to address the discouraging negative result of Long and Servedio
(2010), van Rooyen et al. (2015) considered a weakening of the (Long and
Servedio, 2010) conditions for a convex potential function. In particular,
they allow such functions $\phi$ to take negative values (which is disallowed
by the definition of Long and Servedio). We refer to a function satisfying the
condition of van Rooyen et al. (2015) as a “relaxed convex potential
function.”

The main result of van Rooyen et al. (2015) is that they propose a certain
relaxed convex potential function, which we denote $\phi^{\star}$, and prove
that it is “RCN-robust”.111van Rooyen et al. (2015) uses the term “SLN-
robust”, where the acronym stands for Symmetric Label Noise. We give a formal
definition of RCN-robustness in Section 2, but intuitively it means that a
minimizer of this potential function (minimizing over all linear combinations
of base features) under random classification noise performs no worse than a
minimizer obtained with no random classification noise. van Rooyen et al.
(2015) also define a notion of “strong RCN-robustness” and show that their
$\phi^{\star}$ is the unique relaxed convex potential function which satisfies
strong RCN-robustness.

This note. The purpose of the present note is to discuss the _accuracy_ of the
classifier obtained by minimizing the relaxed convex potential function
$\phi^{\star}$ of van Rooyen et al. (2015). Our main observation is that, for
a simple learning problem corresponding to a certain “clean” data distribution
${\cal D}_{2}$ that is linearly separable with a margin, minimizing
$\phi^{\star}$ over all bounded-norm linear combinations of base features
_even when there is no random classification noise_ only yields a binary
classifier with an error rate of 1/2. Since, as shown by van Rooyen et al.
(2015), $\phi^{\star}$ is the unique strong RCN-robust relaxed convex
potential function, this means that minimizing any strong RCN-robust relaxed
convex potential function in this noise-free scenario may only yield a binary
classifier with an error rate of 1/2, which can be obtained through random
guessing.

Our observation is consistent with the result of van Rooyen et al. (2015) that
$\phi^{\star}$ is RCN-robust, since, informally, that condition only states
that “you don’t do any worse when there is RCN than when there is no RCN.” Our
example demonstrates even when there is no noise, the accuracy of the binary
classifier obtained by minimizing $\phi^{\star}$ may be only 1/2, and of
course the accuracy is no worse than this when there actually is noise.

##  2 Preliminaries

###  2.1 Background: the negative result of Long and Servedio (2010) for
convex potential functions

Convex potential functions. We recall the following definition which is
central to the work of Long and Servedio (2010):

######  Definition 1 ((Long and Servedio, 2010), Definition 1).

A function $\phi:\mathds{R}\to\mathds{R}$ is a _convex potential function_ if
it satisfies the following:

  1. 1.

$\phi\in C^{1}$ (i.e. $\phi$ is differentiable and $\phi^{\prime}$ is
continuous) and $\phi$ is convex and nonincreasing; and

  2. 2.

$\phi^{\prime}(0)<0$ and $\lim_{x\to\infty}\phi(x)=0$ (hence $\phi$ is
everywhere non-negative).

A number of potential functions used in the literature fit this definition,
including the exponential potential function used by AdaBoost (Freund and
Schapire, 1997), the mixed linear/exponential potential function used by
MadaBoost (Domingo and Watanabe, 2000), and the logistic function used by
LogitBoost (Friedman et al., 1998); see Table 1.

Potential function |  Reference |  Satisfies Definition 1?
---|---|---
Exponential: |  |
$\phi(z)=e^{-z}$ |  (Freund and Schapire, 1997) |  Yes
Mixed linear/exponential: |  |
$\phi(z)=\begin{cases}1-z&\text{~{}if~{}}z\leq 0\\\ e^{-z}&\text{~{}if~{}}z>0\\\ \end{cases}$ |  (Domingo and Watanabe, 2000) |  Yes
Logistic: |  |
$\phi(z)=\ln(1+e^{-2z})$ |  (Friedman et al., 1998) |  Yes
Hinge: |  |
$\phi(z)=\max\\{0,1-z\\}$ |  (Gentile and Warmuth, 1998) |  No
Unhinged: |  |
$\phi(z)=1-z$ |  (van Rooyen et al., 2015) |  No
Table 1: Some commonly used potential functions.

Linearly separable learning problems. One of the simplest models for binary-
labeled data over $\mathds{R}^{d}$ is that of data which is _linearly
separable with a margin_. A “clean” probability distribution ${\cal D}$ over
$\mathds{R}^{d}\times\\{-1,1\\}$ is linearly separable with margin $\gamma>0$
if there is a target weight vector $w=(w_{1},\dots,w_{d})\in\mathds{R}^{d}$
such that

$\mathop{{\bf Pr}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}\left[{\frac{\bm{y}(w\cdot\bm{x})}{|w_{1}|+\cdots+|w_{d}|}}<\gamma\right]=0.$

A very standard learning approach for such a setting is to choose a hypothesis
weight vector $v=(v_{1},\dots,v_{d})\in\mathds{R}^{d}$ with the aim of
minimizing the “global” potential function

$P_{\phi,{\cal D}}(v):=\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right].$ (1)

(Of course, given a finite sample of draws from ${\cal D}$, this is typically
done by minimizing the corresponding expectation over the sample.) It is well
known that if ${\cal D}$ is linearly separable with margin $\gamma>0$, then
for a range of different choices of the convex potential function $\phi$
(including the AdaBoost, MadaBoost and LogitBoost potential functions
described above), greedy iterative algorithms that perform coordinatewise
gradient descent to minimize $P_{\phi,{\cal D}}$ will drive the
misclassification error $\operatorname{{\bf Pr}}_{(\bm{x},\bm{y})\sim{\cal
D}}[\bm{y}\neq\mathrm{sign}(v\cdot\bm{x})]$ to zero. Indeed, the AdaBoost
(Freund and Schapire, 1997), MadaBoost (Domingo and Watanabe, 2000) and
LogitBoost (Friedman et al., 1998) boosting algorithms correspond precisely to
greedy coordinatewise gradient descent procedures of this sort; see the work
of Mason et al. (1999) for details.

Learning problems with random classification noise. Let ${\cal D}$ be a data
distribution over $\mathds{R}^{d}\times\\{-1,1\\}$ as described above. The
_$\eta$ -RCN corrupted_ version of ${\cal D}$ is the following distribution
$\overline{{\cal D}}_{\eta}$ over $\mathds{R}^{d}\times\\{-1,1\\}$: a draw
from $\overline{{\cal D}}_{\eta}$ is obtained by drawing
$(\bm{x},\bm{y})\sim{\cal D}$ and flipping the label $\bm{y}$ with probability
$\eta$.

The negative result of Long and Servedio (2010). The main result of Long and
Servedio (2010) is that there is _no_ convex potential function such that
minimizing $\phi$ on the $\eta$-RCN corrupted distribution $\overline{\cal
D}_{\eta}$ will succeed in achieving nontrivial classification accuracy:

######  Theorem 2.

Fix any noise rate $0<\eta<1/2$ and any convex potential function $\phi$.
There is a distribution ${\cal D}$ over $\mathds{R}^{2}\times\\{-1,1\\}$ (in
fact the distribution ${\cal D}$ is supported on three points in the unit
disc) and a margin parameter $\gamma>0$ such that (a) ${\cal D}$ is linearly
separable with margin $\gamma$, but (b) any weight vector $v$ which minimizes
$\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim\overline{{\cal
D}}_{\eta}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$ has $\operatorname{{\bf
Pr}}_{(\bm{x},\bm{y})\sim{\cal
D}}[\bm{y}\neq\mathrm{sign}(v\cdot\bm{x})]=1/2.$

(See Appendix A for a proof that for any convex potential function $\phi$, the
minimizer analyzed in Theorem 2 exists.)

###  2.2 Relaxed convex potential functions: a new hope?

Motivated by the goal of circumventing the negative result of Theorem 2, van
Rooyen et al. (2015) consider a relaxed form of Definition 1:

######  Definition 3.

A function $\phi:\mathds{R}\to\mathds{R}$ is a _relaxed_ convex potential
function if it satisfies the following:

  1. 1.

$\phi\in C^{1}$ (i.e. $\phi$ is differentiable and $\phi^{\prime}$ is
continuous) and $\phi$ is convex and nonincreasing; and

  2. 2.

$\phi^{\prime}(0)<0$.

The only difference between Definition 1 and Definition 3 is that the latter
does not require $\lim_{x\to\infty}\phi(x)=0$; a relaxed convex loss function
may take (arbitrarily large magnitude) negative values. van Rooyen et al.
(2015) exploit this flexibility by proposing the following simple potential
function, which they call the “unhinged loss”:

$\phi^{\ast}(z)=1-z.$

It is trivial to verify that $\phi^{\star}$ satisfies Definition 3 and hence
is a valid relaxed convex potential function. (Note, also, that if the simpler
$\phi(z)=-z$ is used instead, all gradients and minima are unaffected.) We
note that the unhinged loss is a member of the class of _symmetric_ loss
functions, which satisfy $\phi(z)+\phi(-z)=$constant; such loss functions have
been studied by a number of authors, see e.g. Charoenphakdee et al. (2019a);
Ghosh et al. (2015).

RCN-robustness. van Rooyen et al. (2015) analyze the relaxed convex potential
function $\phi^{\star}$ through the lens of a new notion which we will call
_RCN-robustness_. Their definition (Definition 1 of van Rooyen et al. (2015))
applies to a general pair $(\ell,{\cal F})$ where $\ell$ is a loss function
and ${\cal F}$ is a class which may consist of any collection of functions
mapping a domain $X$ to $\mathds{R}$.

Informally, a pair $(\phi,{\cal F})$ is RCN-robust if minimizing $\phi$ over
${\cal F}$ on noise-free data gives the same binary classification performance
as minimizing $\phi$ over ${\cal F}$ on RCN-contaminated data. More precisely,
we have the following:

######  Definition 4.

Let ${\cal F}$ be a set of real-valued functions over $\mathds{R}^{d}$ and let
$\phi$ be a potential function. The pair $(\phi,{\cal F})$ is is said to be
_RCN-robust_ if the following holds: Let ${\cal D}$ be any distribution over
$\mathds{R}^{d}\times\\{-1,1\\}$ and let $0<\eta<1/2$ be any noise rate. If
$f$ is a minimizer of $\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}\left[\phi(\bm{y}(f(\bm{x})))\right]$ over $f\in{\cal F}$ and $g$ is the
minimizer of $\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim\overline{\cal
D}_{\eta}}\left[\phi(\bm{y}(g(\bm{x})))\right]$ over $g\in{\cal F}$, then

$\mathop{{\bf Pr}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}[\bm{y}\neq\mathrm{sign}(f(\bm{x}))]=\mathop{{\bf
Pr}\/}_{(\bm{x},\bm{y})\sim{\cal D}}[\bm{y}\neq\mathrm{sign}(g(\bm{x}))].$ (2)

van Rooyen et al. (2015) specialize Definition 4 to the function class ${\cal
F}_{\mathrm{lin}}$ of all linear functions $x\mapsto v\cdot x$ from
$\mathds{R}^{d}\to\mathds{R}$ (see Section 3.2 of their paper). However, a
problem with Definition 4 for this function class is that $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$
may not have a minimum. This is not merely a technicality. In fact, for
standard loss functions such as the logistic loss or the exponential loss, for
any linearly separable distribution $\cal D$, $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$
does not have a minimum, informally, because scaling up $v$ increases all of
the margins, which decreases all of the losses.222Implicit bias research
analyzes the effect of the algorithm that drives $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$ to
zero on the classification behavior of the limiting classifier. Different
algorithms lead to markedly different limiting classifiers (Telgarsky, 2013;
Soudry et al., 2018; Ji and Telgarsky, 2019) .  van Rooyen et al. (2015)
interpret Theorem 2 as saying that for $d\geq 2$, the pair $(\phi,{\cal
F}_{\mathrm{lin}})$ cannot be RCN-robust for any convex potential function
$\phi$ (see Proposition 1 of Section 3.2 of their paper), but the fact that
the minimizer typically doesn’t exist in the absence of noise interferes with
this interpretation. The unhinged loss also cannot be minimized over ${\cal
F}_{\mathrm{lin}}$, since by scaling up the weight vector of any linear
separator, the unhinged loss can achieve an arbitrarily large negative value.

van Rooyen et al. (2015) also consider the class ${\cal F}_{\mathrm{lin},r}$
of all linear functions whose weight vector has length at most $r$. They prove
the following:

######  Theorem 5 ((van Rooyen et al., 2015), Section 5.1).

For all $d$ and all $r>0$, $(\phi^{\star},{\cal F}_{\mathrm{lin},r})$ is RCN-
robust.

van Rooyen et al. (2015) further establish a number of additional properties
about the unhinged loss $\phi^{\star}$; most of these will not concern us, but
one simple property, which we now explain, is relevant to our discussion in
Section 3. As above let $v\in\mathds{R}^{d}$ be the minimizer of $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$
subject to $||v||\leq r$ and let $v^{\prime}\in\mathds{R}^{d}$ be the
minimizer of $\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim\overline{\cal
D}_{\eta}}\left[\phi(\bm{y}(v^{\prime}\cdot\bm{x}))\right]$ subject to
$||v^{\prime}||\leq r$. van Rooyen et al. (2015) make the straightforward but
useful observation that $v$ is the vector corresponding to a “nearest centroid
classifier” (see Servedio (2002), (Tibshirani et al., 2002), p. 181 of
(Manning et al., 2008), and Section 5.1 of Shawe-Taylor and Cristianini
(2004)), i.e. we have

$v=\alpha\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}[\bm{y}\bm{x}]$ (3)

for a suitable rescaling factor $\alpha$, and furthermore that $v=v^{\prime}$
(this holds since the values of $\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}\left[\phi^{\star}(\bm{y}(w\cdot\bm{x}))\right]$ and $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim\overline{{\cal
D}}_{\eta}}\left[\phi^{\star}(\bm{y}(w\cdot\bm{x}))\right]$ are linearly
related with a slope of $1-2\eta>0$).

##  3 A separable learning problem where minimizing the unhinged loss on clean
data yields a poor classifier

$v$(1,0)$(\gamma,\sqrt{1-\gamma^{2}})$$(\gamma,-2\gamma)$ Figure 1: The
distribution ${\cal D}$ over $\mathds{R}^{2}\times\\{-1,1\\}$. All examples
have label $+1$; the two examples with weight $1/4$ are depicted with small
filled circles and the example with weight $1/2$ is depicted with a larger
unfilled circle.

In this section we observe that while the unhinged loss $\phi^{\star}$ is
strongly robust, there are simple linearly separable data distributions for
which minimizing $\phi^{\star}$ over all functions in ${\cal
F}_{\mathrm{lin},r}$ _even in the absence of random classification noise_ only
yields a binary classifier with an error rate of 1/2. So while (2) is
satisfied, it holds because both error rates are equal to 1/2.

We illustrate this with the linearly separable learning scenario which is
depicted in Figure 1. The distribution ${\cal D}$ over
$\mathds{R}^{2}\times\\{-1,1\\}$ is as follows: given a parameter
$0<\gamma<0.0901$,

  * •

${\cal D}$ puts weight $1/4$ on the labeled example
$x^{(1)}:=(1,0),y^{(1)}:=1$;

  * •

${\cal D}$ puts weight $1/4$ on the labeled example
$x^{(2)}:=(\gamma,\sqrt{1-\gamma^{2}}),y^{(2)}:=1$;

  * •

${\cal D}$ puts weight $1/2$ on the labeled example
$x^{(3)}:=(\gamma,-2\gamma),y^{(3)}:=1$.

It is clear that ${\cal D}$ is linearly separable with margin $\gamma$. By
(3), the vector in $\mathds{R}^{2}$ which minimizes $\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right]$
points in the direction of

$v=(v_{1},v_{2})=\left({\frac{1}{4}}+{\frac{3\gamma}{4}},{\frac{\sqrt{1-\gamma^{2}}}{4}}-\gamma\right).$

For $0<\gamma<0.0901$ we have $v\cdot x^{(3)}<0$ and hence
$\mathrm{sign}(v\cdot x^{(3)})\neq y^{(3)}$, so the LHS of (2) is $1/2$. Since
$v^{\prime}=v$ the RHS of (2) is also $1/2.$

##  4 Implicit bias

This section includes a couple of observations about the implicit bias of
algorithms that iteratively reduce the unhinged loss. Analogous results have
been obtained for other loss functions (Telgarsky, 2013; Soudry et al., 2018;
Ji and Telgarsky, 2019).

###  4.1 Gradient descent

Recall that the unhinged loss function is defined to be $\phi^{\ast}(z)=1-z.$
If ${\cal D}$ is uniform over
$(x_{1},y_{1}),...,(x_{n},y_{n})\in\mathds{R}^{d}\times\\{-1,1\\}$, then for
any $v\in\mathds{R}^{d}$ the gradient of the unhinged loss at $v$ is
$-\sum_{i=1}^{n}y_{i}x_{i}$ (note that this does not depend on $v$). Thus, if
the unhinged loss is minimized by gradient descent starting with an initial
solution of $0$, all iterates are multiples of $\sum_{i=1}^{n}y_{i}x_{i}$. If
the initial solution is $v_{0}$, then, after $T$ updates with step size
$\eta$, the weight vector is $v_{0}+\eta T\sum_{i=1}^{n}y_{i}x_{i}$. As $T$
goes to infinity, the angle between this weight vector and
$\sum_{i=1}^{n}y_{i}x_{i}$ goes to zero.

###  4.2 Coordinate descent

As mentioned earlier, popular boosting algorithms can be viewed as coordinate
descent on a convex potential function, which works by repeatedly finding the
coordinate axis with the steepest descent direction and making an update in
that direction. Informally, the unhinged loss rewards increasing the margin
$yv\cdot x$ of a correctly classified example $(x,y)$ as much as increasing
the negative margin of an incorrectly classified example, but increasing the
margin of a correctly classified example does not make progress towards
overall classification accuracy. (In contrast, the tendency of the exponential
loss to place more importance on misclassified examples is key to AdaBoost’s
ability to boost.) If we denote the components of $x_{i}$ by
$x_{i,1},...,x_{i,d}$, since the gradient of the unhinged loss is the same for
all $v$, when it is minimized by coordinate descent starting from the zero
weight vector all of its iterates will only have nonzero components on members
of $\mathrm{argmax}_{j}\sum_{i}y_{i}x_{ij}$. (If there is not a tie for the
best weak learner, this will be a single component.) From a boosting point of
view, a boosting algorithm based on the unhinged loss allows a weak learner to
keep returning the same (weak) hypothesis.

##  5 Discussion

While the unhinged loss is noise-tolerant in a sense, minimizing it can fail
to find an accurate classifier on data that is linearly separable with a large
margin. On the other hand, minimizing the unhinged loss has been found to
yield reasonable accuracy on natural data (see Patrini et al., 2017;
Charoenphakdee et al., 2019b). This is not entirely unexpected, since, when
learning linear models, minimizing the unhinged loss is closely related to
performing Naive Bayes classification (Domingos and Pazzani, 1997; Ng and
Jordan, 2001), using a spherical Gaussian to model the class-conditional
distributions.

Given our results, one natural goal for future work is to study whether there
are conditions on potential functions which achieve an attractive tradeoff
between noise-robustness and usefulness for learning (in the sense that
minimizing the potential function yields an accurate classifier). Tools
developed for studying Fisher consistency (Fisher, 1922), consistent loss
functions (Zhang, 2004), classification calibration (Bartlett et al., 2006)
and $H$-consistency (Long and Servedio, 2013) may be useful for this. In
particular, it would be interesting to investigate symmetric potential
functions (see e.g. Charoenphakdee et al. (2019a); Ghosh et al. (2015)) and
the multiclass setting (see e.g. (Ghosh et al., 2017; Zhang and Sabuncu,
2018)) in light of this question.

##  Appendix A A minimizer exists for noisy data

In this appendix we show that for all convex potential functions $\phi$, all
finite-covariance distributions ${\cal D}$ over $\mathds{R}^{d}$, and all
$\eta\in(0,1/2)$, the function $P_{\phi,\overline{{\cal D}}_{\eta}}$ has a
minimum. (Recall from (1) that $P_{\phi,\overline{{\cal
D}}_{\eta}}(v):=\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim\overline{{\cal
D}}_{\eta}}\left[\phi(\bm{y}(v\cdot\bm{x}))\right].$)

We recall some useful background.

######  Definition 6.

For any $a$, the set $\\{x:f(x)\leq a\\}$ is a level set for
$f:\mathds{R}^{d}\to\mathds{R}$.

######  Definition 7.

A nonzero vector $u\in\mathds{R}^{d}$ is a direction of recession for a
function $f$ if, for all nonempty level sets $L$ of $f$, there exists some
$x_{0}$ such that $x_{0}+\lambda u\in L$ for all $\lambda\geq 0$. (Informally,
all non-empty level sets of $f$ extend infinitely in the $u$ direction.)

######  Lemma 8 ((Rockafellar, 2015), Theorem 27.1).

The set of minima of a continuous convex function $f$ is nonempty and bounded
iff $f$ does not have any direction of recession.

######  Definition 9.

Say that a probability distribution ${\cal D}$ over
$\mathds{R}^{d}\times\\{-1,1\\}$ has finite covariance if, for all unit length
$u\in\mathds{R}^{d}$, $\operatorname{{\bf E}}_{(\bm{x},\bm{y})\sim{\cal
D}}[(u\cdot\bm{x})^{2}]$ exists.

Now we are ready to analyze $P_{\phi,\overline{{\cal D}}_{\eta}}$.

######  Proposition 10.

For any convex potential function $\phi$, for any $\eta\in(0,1/2)$, for any
finite-covariance distribution ${\cal D}$ over $\mathds{R}^{d}$,
$P_{\phi,\overline{{\cal D}}_{\eta}}$ has a minimum.

###### Proof.

First, we may assume without loss of generality that, for all unit length
$u\in\mathds{R}^{d}$,

$\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}[|u\cdot\bm{x}|]>0$ (4)

since otherwise $u\cdot\bm{x}=0$ almost surely, and $P_{\phi,\overline{{\cal
D}}_{\eta}}(v)$ is unaffected by projecting $v$ onto the subspace of
$\mathds{R}^{d}$ orthogonal to $u$.

Assume for contradiction that some $u$ is a direction of recession for
$P_{\phi,\overline{{\cal D}}_{\eta}}$. For any $x_{0}\in\mathds{R}^{d}$ and
$\lambda\geq 0$, we have

$\displaystyle P_{\phi,\overline{{\cal D}}_{\eta}}(x_{0}+\lambda u)$
$\displaystyle=\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim\overline{{\cal
D}}_{\eta}}[\phi(\bm{y}((x_{0}+\lambda u)\cdot\bm{x}))]$ (def. of
$P_{\phi,\overline{{\cal D}}_{\eta}}$) $\displaystyle=\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal D}}[\eta\phi(-\bm{y}((x_{0}+\lambda
u)\cdot\bm{x}))+(1-\eta)\phi(\bm{y}((x_{0}+\lambda u)\cdot\bm{x}))]$ (def. of
$\overline{{\cal D}}_{\eta}$) $\displaystyle\geq\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}[\eta\max_{\tilde{y}\in\\{-1,1\\}}\phi(\tilde{y}((x_{0}+\lambda
u)\cdot\bm{x}))]$ (since $\eta<1/2$ and $\phi\geq 0$)
$\displaystyle=\mathop{{\bf E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}[\eta\phi(-|(x_{0}+\lambda u)\cdot\bm{x}|)]$ (since $\phi$ is
nonincreasing) $\displaystyle\geq\eta\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}[\phi(0)-\phi^{\prime}(0)\cdot|(x_{0}+\lambda u)\cdot\bm{x}|]$ (since
$\phi$ is convex and $\phi^{\prime}(0)<0$)
$\displaystyle\geq\eta\left(\phi(0)-\phi^{\prime}(0)\lambda\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})\sim{\cal
D}}[|u\cdot\bm{x}|]+\phi^{\prime}(0)\mathop{{\bf
E}\/}_{(\bm{x},\bm{y})}[|x_{0}\cdot\bm{x}|]\right).$ (triangle inequality,
$\phi^{\prime}(0)<0$)

Thus $\lim_{\lambda\rightarrow\infty}P_{\phi,\overline{{\cal
D}}_{\eta}}(x_{0}+\lambda u)=\infty$, which contradicts the assumption that
$u$ is a direction of recession for $P_{\phi,\overline{{\cal D}}_{\eta}}$. ∎

## References

  * Bartlett et al. [2006] Peter L Bartlett, Michael I Jordan, and Jon D McAuliffe.  Convexity, classification, and risk bounds.  _Journal of the American Statistical Association_ , 101(473):138–156, 2006.
  * Charoenphakdee et al. [2019a] Nontawat Charoenphakdee, Jongyeong Lee, and Masashi Sugiyama.  On symmetric losses for learning from corrupted labels.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 961–970. PMLR, 2019a.
  * Charoenphakdee et al. [2019b] Nontawat Charoenphakdee, Jongyeong Lee, and Masashi Sugiyama.  On symmetric losses for learning from corrupted labels.  In _International Conference on Machine Learning_ , pages 961–970. PMLR, 2019b.
  * Dietterich [2000] T.G. Dietterich.  An experimental comparison of three methods for constructing ensembles of decision trees: bagging, boosting, and randomization.  _Machine Learning_ , 40(2):139–158, 2000.
  * Domingo and Watanabe [2000] C. Domingo and O. Watanabe.  MadaBoost: a modified version of AdaBoost.  In _Proceedings of the Thirteenth Annual Conference on Computational Learning Theory (COLT)_ , pages 180–189, 2000.
  * Domingos and Pazzani [1997] Pedro Domingos and Michael Pazzani.  On the optimality of the simple bayesian classifier under zero-one loss.  _Machine learning_ , 29(2):103–130, 1997.
  * Fisher [1922] Ronald A Fisher.  On the mathematical foundations of theoretical statistics.  _Philosophical transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character_ , 222(594-604):309–368, 1922.
  * Freund and Schapire [1996] Y. Freund and R. Schapire.  Experiments with a new boosting algorithm.  In _Proceedings of the Thirteenth International Conference on Machine Learning_ , pages 148–156, 1996.
  * Freund and Schapire [1997] Y. Freund and R. Schapire.  A decision-theoretic generalization of on-line learning and an application to boosting.  _Journal of Computer and System Sciences_ , 55(1):119–139, 1997.
  * Friedman et al. [1998] J. Friedman, T. Hastie, and R. Tibshirani.  Additive logistic regression: a statistical view of boosting.  _Annals of Statistics_ , 28(2):337–407, 1998\.
  * Gentile and Warmuth [1998] Claudio Gentile and Manfred K. Warmuth.  Linear hinge loss and average margin.  In _Advances in Neural Information Processing Systems 11_ , pages 225–231. The MIT Press, 1998.
  * Ghosh et al. [2015] Aritra Ghosh, Naresh Manwani, and P. S. Sastry.  Making risk minimization tolerant to label noise.  _Neurocomputing_ , 160:93–107, 2015.
  * Ghosh et al. [2017] Aritra Ghosh, Himanshu Kumar, and PS Sastry.  Robust loss functions under label noise for deep neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Ji and Telgarsky [2019] Ziwei Ji and Matus Telgarsky.  The implicit bias of gradient descent on nonseparable data.  In _Conference on Learning Theory_ , pages 1772–1798. PMLR, 2019\.
  * Long and Servedio [2013] Phil Long and Rocco Servedio.  Consistency versus realizable h-consistency for multiclass classification.  In _International Conference on Machine Learning_ , pages 801–809. PMLR, 2013.
  * Long and Servedio [2010] Philip M Long and Rocco A Servedio.  Random classification noise defeats all convex potential boosters.  _Machine learning_ , 78(3):287–304, 2010.
  * Maclin and Opitz [1997] R. Maclin and D. Opitz.  An empirical evaluation of bagging and boosting.  In _Proceedings of the Fourteenth National Conference on Artificial Intelligence and Ninth Innovative Applications of Artificial Intelligence Conference (AAAI/IAAI)_ , pages 546–551, 1997.
  * Manning et al. [2008] Christopher D. Manning, Prabhakar Raghavan, and Hinrich Schütze.  _Introduction to Information Retrieval_.  Cambridge University Press, 2008.
  * Manwani and Sastry [2013] Naresh Manwani and P. S. Sastry.  Noise tolerance under risk minimization.  _IEEE Trans. Cybern._ , 43(3):1146–1151, 2013\.
  * Mason et al. [1999] L. Mason, J. Baxter, P. Bartlett, and M. Frean.  Boosting algorithms as gradient descent.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 512–518, 1999.
  * Natarajan et al. [2013] Nagarajan Natarajan, Inderjit S. Dhillon, Pradeep Ravikumar, and Ambuj Tewari.  Learning with noisy labels.  In _Advances in Neural Information Processing Systems 26: 27th Annual Conference on Neural Information Processing Systems 2013_ , pages 1196–1204, 2013.
  * Ng and Jordan [2001] A. Y. Ng and M. I. Jordan.  On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes.  _NIPS_ , 2001.
  * Patrini et al. [2017] Giorgio Patrini, Alessandro Rozza, Aditya Krishna Menon, Richard Nock, and Lizhen Qu.  Making deep neural networks robust to label noise: A loss correction approach.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1944–1952, 2017.
  * Rockafellar [2015] Ralph Tyrell Rockafellar.  _Convex analysis_.  Princeton university press, 2015.
  * Servedio [2002] Rocco A. Servedio.  Perceptron, winnow, and PAC learning.  _SIAM J. Comput._ , 31(5):1358–1369, 2002.
  * Shawe-Taylor and Cristianini [2004] Jonathan Shawe-Taylor and Nello Cristianini.  _Kernel Methods for Pattern Analysis_.  Cambridge University Press, 2004.
  * Soudry et al. [2018] Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  _The Journal of Machine Learning Research_ , 19(1):2822–2878, 2018.
  * Telgarsky [2013] Matus Telgarsky.  Margins, shrinkage, and boosting.  In _ICML_ , volume 28 of _JMLR Workshop and Conference Proceedings_ , pages 307–315, 2013.
  * Tibshirani et al. [2002] Robert Tibshirani, Trevor Hastie, Balasubramanian Narasimhan, and Gilbert Chu.  Diagnosis of multiple cancer types by shrunken centroids of gene expression.  _Proceedings of the National Academy of Sciences_ , 99(10):6567–6572, 2002.
  * van Rooyen et al. [2015] Brendan van Rooyen, Aditya Menon, and Robert C. Williamson.  Learning with Symmetric Label Noise: The Importance of Being Unhinged.  In _Advances in Neural Information Processing Systems 28_ , 2015.
  * Zhang [2004] Tong Zhang.  Statistical behavior and consistency of classification methods based on convex risk minimization.  _The Annals of Statistics_ , 32(1):56–85, 2004\.
  * Zhang and Sabuncu [2018] Zhilu Zhang and Mert Sabuncu.  Generalized cross entropy loss for training deep neural networks with noisy labels.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.


††thanks: $\dagger$NASA Hubble Fellow

# Does God play dice with star clusters?

Michael Y. Grudić1,† Stella S. R. Offner2 Dávid Guszejnov3,† Claude-André
Faucher-Giguère4 Philip F. Hopkins5 1Carnegie Observatories, 813 Santa Barbara
St, Pasadena, CA 91101, USA  2Department of Astronomy, The University of Texas
at Austin, TX 78712, USA  3Center for Astrophysics $|$ Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, USA  4CIERA and Department of Physics
and Astronomy, Northwestern University, 1800 Sherman Ave, Evanston, IL 60201,
USA  5TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA  Corresponding author<EMAIL_ADDRESS>
###### Abstract

When a detailed model of a stellar population is unavailable, it is most
common to assume that stellar masses are independently and identically
distributed according to some distribution: the universal initial mass
function (IMF). However, stellar masses resulting from causal, long-ranged
physics cannot be truly random and independent, and the IMF may vary with
environment. To compare stochastic sampling with a physical model, we run a
suite of 100 STARFORGE radiation magnetohydrodynamics simulations of low-mass
star cluster formation in $2000M_{\odot}$ clouds that form $\sim 200$ stars
each on average. The stacked IMF from the simulated clouds has a sharp
truncation at $\sim 28M_{\odot}$, well below the typically-assumed maximum
stellar mass $M_{\rm up}\sim 100-150M_{\odot}$ and the total cluster mass. The
sequence of star formation is not totally random: massive stars tend to start
accreting sooner and finish later than the average star. However, final
cluster properties such as maximum stellar mass and total luminosity have a
similar amount of cloud-to-cloud scatter to random sampling. Therefore
stochastic sampling does not generally model the stellar demographics of a
star cluster as it is forming, but may describe the end result fairly well, if
the correct IMF – and its environment-dependent upper cutoff – are known.

###### keywords:

stars: mass function – galaxies: star clusters: general – stars: formation

##  1 Introduction

Astrophysicists aim to interpret light from galaxies and star clusters in
terms of their stellar populations and the physical processes that affect
them. More often than not, detailed knowledge of the individual stars is not
available, so inferences about the system rely on a model of their stellar
demographics. The cornerstone of such models is the stellar initial mass
function (IMF): the distribution of stellar birth masses (Salpeter, 1955;
Hopkins, 2018). The standard approach to modeling a galaxy or star cluster is
to sample stellar masses from the IMF and evolve the stars forward in time
according to a stellar evolution model (e.g. Leitherer et al., 1999; Bruzual &
Charlot, 2003).

A common model assumption is that stellar masses are sampled independently and
identically from an invariant IMF (“stochastic sampling”), i.e., the
probability of a star forming with mass $m_{\rm 1}$ does not depend on the
mass of another star $m_{\rm 2}$, nor the macroscopic properties of the system
the star is forming in. Or, in other words, 10 little star clusters are
statistically identical to one big one. This is the simplest approach and is
motivated by the observed lack of strong, systematic variations in the IMF in
normal star-forming environments (Bastian et al., 2010; Offner et al., 2014).

Stochastic sampling should not necessarily be interpreted physically: the
initial gas flows in molecular clouds that produce individual stars cannot be
truly random and independent of each other. The formation of one star can
affect another via gravity and feedback in the form of protostellar outflows,
winds, radiation, supernovae, and cosmic rays. So stochastic sampling cannot
be the ground-truth star formation process, and if star clusters deviate from
it significantly then so will any resulting inferences from the integrated
light of galaxies (Weidner & Kroupa, 2006; Kroupa et al., 2013; Kroupa &
Jerabkova, 2021).

Of course, all models are wrong at some level, but some are useful. Stochastic
sampling from an invariant IMF has enjoyed wide application due to its
simplicity and the lack of an alternative model with compelling physical
motivation. It has been used to characterize star clusters and their
demographics from photometry (Fumagalli et al., 2011; da Silva et al., 2012),
where it gives a marginally-better fit to the data than IMF-averaged
photometric models without stochasticity (Krumholz et al., 2015) and can
enable inferences about the highly-incomplete low-mass cluster population, if
accurate (Tang et al., 2023). Stochastic sampling has also seen increasing use
in high-resolution galaxy and molecular cloud simulations (Fujii & Portegies
Zwart, 2015; Hu et al., 2017; Su et al., 2018; Lahén et al., 2020; Smith,
2021; Liow et al., 2022). On molecular cloud ($\lesssim 100\rm pc$) scales,
the granularity of feedback and gravity from individual stars is quite
important; the specific timing, sampling, and distribution of massive star
formation has a large effect on subsequent star formation (Grudić & Hopkins,
2019; Dinnbier & Walch, 2020; Lewis et al., 2023). However, lacking a detailed
physical picture of star cluster formation, it has been difficult to assess
the fidelity of stochastic sampling fully.

Recently, the STARFORGE framework111https://www.starforge.space (Grudić et
al., 2021) has enabled radiation magnetohydrodynamics (RMHD) star formation
simulations with significantly enhanced realism and self-consistency. By
accounting for key feedback processes and resolving individual star formation,
the simulations can follow the entire sequence of star formation from initial
cloud collapse, through star formation and accretion, feedback, to the
eventual cloud disruption (Grudić et al., 2022). This allows us to predict the
final masses of stars and hence the IMF in a consistent framework (Guszejnov
et al., 2022).

In this paper we present the results of a statistical ensemble of 100
STARFORGE simulations and measure their combined and individual stellar mass
statistics. The statistical power of this suite allows us to compare
stochastic sampling models to a physical model of star cluster formation with
all relevant processes for the first time. In particular we aim to test
whether statistical sampling, independent of cloud mass, provides a good
representation of star formation process realized in nature. We will focus on
the statistics of massive ($\gtrsim 8M_{\odot}$) stars, which dominate the
light and stellar feedback from young star clusters, and hence are most
relevant for photometric modeling and numerical simulations.

##  2 Simulations

We use the GIZMO code (Hopkins, 2015) to run an ensemble of 100 STARFORGE RMHD
star cluster formation simulations. The simulations use the “full” STARFORGE
physics setup introduced in Grudić et al. (2022), with numerical methods for
sink particles, stellar evolution, and feedback coupling described in full in
Grudić et al. (2021). The simulations solve the equations of ideal
magnetohydrodynamics (MHD) using the GIZMO Meshless Finite Mass MHD solver
(Hopkins & Raives, 2016) with constrained-gradient divergence minimization
(Hopkins, 2016).

Gas cells turn into sink particles if they satisfy all of a list of criteria
intended to identify centers of runaway gas collapse. They then accrete other
gas cells that have sufficient density, proximity to the sink, low angular
momentum, and gravitational boundedness (e.g. Bate et al., 1995). Each sink
particle hosts a star that accretes the sink’s gas reservoir smoothly,
evolving in radius and deuterium abundance according to the Offner et al.
(2009) protostellar evolution model. Stars on the main sequence may still
accrete and move along it, and their bolometric luminosity and radius are
given by the Tout et al. (1996) fits. Stars’ emergent spectra are approximated
as a black-body with $T_{\rm eff}=5780\rm
K\left(L_{\star}/R_{\star}^{2}\right)^{1/4}$, and their metallicity- and
luminosity-dependent mass-loss rate and wind velocity are given by the fitting
formulae in Grudić et al. (2021); Grudić et al. (2022). Massive stars can
progress to a Wolf-Rayet phase with strong winds and an eventual supernova,
but this does not occur during the star-forming phase of any of the clouds
simulated here, which are disrupted primarily by protostellar outflows and
ionizing radiation.

We solve the time-dependent, frequency-integrated radiative transfer (RT)
equation with the GIZMO meshless frequency-integrated M1 solver (Hopkins &
Grudić, 2019; Hopkins et al., 2020) with a reduced speed of light of $30\,\rm
km\,s^{-1}$, in 5 frequency bands ranging from Lyman continuum to far IR. The
different bands couple to matter through photoionization, photodissociation,
photoelectric heating, and dust absorption, accounting for photon momentum. We
account for radiative cooling and heating from all major molecular, atomic,
nebular, and continuum processes using the cooling module shared with the
FIRE-3 simulations (Hopkins et al., 2023), except that we perform dust
radiative transfer self-consistently with the IR band of the RT solver,
evolving the dust, gas, and radiation temperatures independently.

###  2.1 Initial conditions

All 100 simulations start with a uniform-density molecular cloud of mass
$M=2000M_{\odot}$ and radius $R=3\rm pc$, in a $30\rm pc$ periodic box filled
with an ambient medium $10^{3}$ times less dense. The initial magnetic field
is uniform and normalized to give a mass-to-flux ratio of 4.2 times the
critical value inside the cloud. The simulations differ only in their initial
random velocity fields, which are different statistical realizations of a
Gaussian random field with a $\propto k^{-2}$ power spectrum, all normalized
so that the initial virial parameter is $\alpha_{\rm
turb}=5\sigma^{2}R/(3GM)=2$, with a natural mix of compressive and solenoidal
modes. We use a fixed quasi-Lagrangian mass resolution of $10^{-3}M_{\odot}$
for normal gas cells and $10^{-4}M_{\odot}$ for protostellar jet and stellar
wind cells, hence each cloud initially consists of $2\times 10^{6}$ gas cells.

Figure 1: Histogram of the 19,009 ZAMS stellar masses resulting from the
stacked 100 star cluster formation simulations. For comparison we plot
estimated $\pm\sigma$ probable ranges of the Kroupa (2001) and Chabrier (2005)
Milky Way field IMFs, using parameter uncertainties quoted in Kroupa (2001)
and assuming $\log m_{\rm c}/M_{\odot}=\log 0.2\pm 0.1\rm dex$ and
$\sigma=0.55\pm 0.10\rm dex$ for the Chabrier (2005) form, controlling all
IMFs for the number of $>0.1M_{\odot}$ stars. The dashed line plots the best-
fit $\propto M^{-2}$ power-law, and the shaded region plots the $\pm 3\sigma$
confidence region of a modified Schechter function fit, sampled from the
posterior distribution of the model parameters given the $>1M_{\odot}$ sample.
The two lower histograms show the IMF obtained from one $2\times
10^{4}M_{\odot}$ cloud run with the same physics model, and a run with the
same initial conditions but no feedback.

##  3 Results

All simulations are run for $10\,\rm Myr$ (5 cloud freefall times) and follow
the usual qualitative sequence predicted by molecular cloud simulations with
star formation and feedback, e.g. as documented in detail in Grudić et al.
(2022). Turbulent motions seed dense substructures that go on to collapse into
stars, which then proceed to accrete and disperse the cloud via feedback,
eventually halting star formation. The clouds form a total of 19,009 stars
totaling $1.26\times 10^{4}M_{\odot}$, for an average final cluster mass of
$126M_{\odot}$. Roughly 80% of the stellar systems appear to be expanding
freely at the end of their respective simulations, and the others have a
steady half-mass radius and are likely bound and destined to dissolve over a
longer timescale. For the purposes of this paper we adopt a broad definition
of “cluster” and do not distinguish between expanding associations and bound
clusters.

The spectrum of zero-age main sequence masses of stars formed in all 100
simulations is plotted in Figure 1, with the Kroupa (2001) and Chabrier (2005)
Solar neighborhood field star IMFs for comparison. This IMF is similar to the
one that we have presented for clouds with these parameters and physics in
previous work (Grudić et al., 2022; Guszejnov et al., 2022), with a power-law
slope of $\sim-2$ in the $1-10M_{\odot}$ range and a turnover or break around
$0.2M_{\odot}$. The median stellar mass is $0.22M_{\odot}$ and the log-
dispersion is $\sigma_{\mathrm{log}M}=0.52\,\rm dex$, similar to estimates for
the Solar neighborhood field population (Chabrier, 2005). The high-mass slope
$\alpha=-2$ is shallower than the normally-assumed -2.35, but well within the
variation measured between individual clusters and within the uncertain range
quoted for the Milky Way field population (Kroupa, 2001).

###  3.1 The high-mass IMF truncation

Stacking 100 realizations unveils a new feature of the IMF: a high-mass
truncation. Figure 1 shows that there is an unambiguous high-mass break around
the sample maximum stellar mass of $28\rm M_{\odot}$, for the 100 low-mass
clusters. To assess this statistically, note that if the $\propto M^{-2}$
power-law portion extrapolated to $M_{\rm up}=150M_{\odot}$, we would expect
$\sim$75 stars more massive than $28M_{\odot}$ in our sample. The probability
of obtaining zero is $2\times 10^{-33}$ assuming Poisson statistics. Even a
steeper $\propto M^{-2.35}$ power-law up to $150M_{\odot}$ can be discarded
with $p=5\times 10^{-12}$ by this argument. The maximum stellar mass found in
these simulations is confidently inconsistent with the size-of-sample effect
from sampling a pure power-law extending to the largest-known stellar masses
(e.g. Crowther et al., 2016), a standard model assumption. The sample maximum
is also significantly smaller than even the smallest cluster mass in the
sample ($\sim 80M_{\odot}$), so it is not directly explained by the limited
mass budget.

There is no physical reason to expect a perfectly-sharp truncation, so we
model the high-mass ($>1M_{\odot}$) IMF as a power-law with a smooth
truncation factor. A modified Schechter function with varying sharpness

$\mathrm{d}N/\mathrm{d}M\propto
M^{\alpha}\exp\left(-\left(M/M^{\ast}\right)^{p}\right)$ (1)

provides a better fit than the standard Schechter (1976) form with $p=1$,
e.g., its additional parameter is strongly preferred by the Schwarz (1978)
information criterion ($\Delta\mathrm{BIC}=-17.5$). Assuming flat priors with
$\alpha\in[0,-4]$, $\log M^{\ast}\in[0,4]$, and $p\in[0,10]$ the posterior
distribution gives marginalized parameters $\alpha=-2.00\pm 0.02$,
$M^{\ast}=24\pm 4M_{\odot}$, and $p=7\pm 2$, and we plot the $3\sigma$
confidence interval of the parameter distribution in Figure 1.

This sample IMF is also statistically distinct from those found in various
previous simulations using the same setup with different initial conditions
(Grudić et al., 2022; Guszejnov et al., 2022). Most of those exhibited the
same high-mass slope of -2 but produced more-massive stars in clusters with
fewer members. For example, in Figure 1 we also plot the IMF from a cloud run
with an identical physics setup, surface density, virial parameter, and mass-
to-flux ratio, but 10 times greater mass ($2\times 10^{4}M_{\odot}$). This
single cloud produces about $10\times$ fewer stars than the 100 smaller
clouds. Despite this, the massive cloud produces 3 stars more massive than the
$28M_{\odot}$ maximum of the smaller clouds, with a maximum of $44M_{\odot}$.
This strongly suggests that the upper cutoff predicted by the STARFORGE models
scales in some way on the properties of the host GMC. We do emphasize that our
model space can still account for the existence of $\gtrsim 100M_{\odot}$
stars - e.g. the $10\times$ higher surface-density cloud in Guszejnov et al.
(2022) formed a $107M_{\odot}$ star.

Figure 1 also shows the IMF from a $2000M_{\odot}$ cloud with the same bulk
properties as the others, but evolved without any stellar feedback physics.
Despite forming only 69 stars, this cloud produced 8 stars more massive than
the most massive star in the feedback simulations, with a maximum stellar mass
of $224M_{\odot}$. Clearly then the upper cutoff is highly sensitive to
feedback, so a dependence on cloud properties is feedback inevitably operates
with variable effectiveness depending on the environment, in both the
momentum-conserving and energy-conserving limits (Fall et al., 2010).

Due to this variation in the upper truncation in the IMF, no sampling model
employing the usual constant $\sim 100-150M_{\odot}$ upper cutoff – stochastic
or not – adequately describes the statistics of the simulated star clusters.
Instead, a model in which the truncation depends on the environment is
required.

###  3.2 Sequence of massive star formation versus random sampling

We will now compare the number of massive ($>8M_{\odot}$) stars in the
clusters at a given time with the mean number expected, $\langle N_{\rm
massive}\rangle=f_{\rm massive}N_{\rm cl}$, where $f_{\rm massive}=0.012$ the
overall fraction of $>8M_{\odot}$ stars in the sample IMF. We count massive
stars in two different ways: (1) we include a protostar in the count after its
formation if it will eventually be massive, and (2) we include a star in the
count only after it has reached $8M_{\odot}$.

The results, binned over the full set of simulation snapshots (over all stages
of cluster formation), are summarized in Figure 2. The average number of
eventually-massive stars tends to be $\sim 3\times$ greater than that expected
from random sampling – i.e., massive stars in the simulations tend to start
forming earlier than expected if stars of different masses formed in random
order. The opposite is true for the counts of presently-massive stars: there
tend to be fewer than expected from random sampling as the cluster forms,
while it has $\lesssim 100$ stars. Stars of different masses do not finish
forming (i.e., reach their maximum mass) at random times. A massive star is
never among the first 20 fully-formed stars in any of these simulations. This
is is likely due to the finite formation time required for the massive stars
to accrete their mass, which can be as long as a few Myr (Grudić et al.,
2022), comparable to the the total lifetime of the simulated clouds.

Figure 2: Expected number of massive ($>8M_{\odot}$) stars in a cluster as a
function of the number of stars in a cluster, measured over time in all
simulations. We count both the number of stars that will eventually exceed
$8M_{\odot}$ and the number of stars that currently exceed $8M_{\odot}$ at a
given number of cluster members. Horizontal bars show bin widths, vertical
bars show counting uncertainties, downward arrows indicate bins with no stars.
For comparison we plot the expected number of massive stars when drawing from
the total sample IMF (Fig. 1) at random (dashed).

###  3.3 Final stellar demographics

Figure 3: Relation between number of stars and maximum stellar mass (left) and
bolometric luminosity (right) at the end of star formation in the 100
simulations; each datapoint represents one simulated star cluster. Error bars
plot the binned median and $\pm\sigma$ range of the simulation results
(horizontal bars show bin widths). For comparison we plot the $\pm\sigma$
quantiles of clusters randomly sampled from sample IMF, a Kroupa (2002) with
standard slope $\alpha=-2.3$ and cutoff $M_{\rm up}=150M_{\odot}$, and a
modified Kroupa IMF with slope $\alpha=-2$ and hard cutoff $M_{\rm
up}=28M_{\odot}$.

Next we examine the statistics of the most important properties for the fully-
formed cluster’s protometric and feedback properties: the maximum stellar mass
$M_{\rm max}$ and the cluster luminosity $L_{\rm bol}$. In Figure 3 we plot
the relation between the number of member stars and $L_{\rm bol}$ and $M_{\rm
max}$ respectively for the simulated clusters. Note that the two plots are
qualitatively identical because $L_{\rm bol}$ is strongly correlated with
$M_{\rm max}$ due to the steep scaling of stellar luminosity with stellar
mass.

The $L_{\rm bol}$ and $M_{\rm max}$ relations exhibit $\sim 0.4\rm dex$ and
$\sim 0.15\rm dex$ of scatter, respectively; the number of stars in the
cluster does not uniquely predict either quantity. Remarkably, the degree of
scatter agrees well with that predicted from stochastic sampling from the
sample IMF, or from a modified Kroupa (2002) IMF with slope $\alpha=-2$ and
hard cutoff $M_{\rm up}=28M_{\odot}$. However, a standard IMF with
$\alpha=-2.3$ and $M_{\rm up}=150M_{\odot}$ exhibits significantly more
scatter in both directions due to the different shape and cutoff.

It is worth re-emphasizing why this result is not necessarily expected: a
priori, the data plotted in Figure 3 did not have to exhibit any significant
scatter whatsoever, let alone resemble random sampling. In principle star
formation could have proceeded as a completely well-determined, organized
process in which the macroscopic system properties tightly constrain the
outcome (e.g. Weidner & Kroupa, 2006). The highly-nonlinear, dissipative
physics at play might have damped out any particulars of the initial seed and
regressed toward an emergent attractor behavior, perhaps a detailed self-
regulation by feedback. But the simulations indicate that this is not the
outcome realized by the physics of star formation in practice.

##  4 Discussion

###  4.1 The variable high-mass truncation of IMF

What does it mean that the most-massive known stars in the Local Group tend to
be found in the most-massive young star clusters? Is it a size-of-sample
effect from sampling a universal power-law, or do certain clusters form under
special conditions that promote the formation of very massive stars? It has
long been known that the relation between cluster mass and maximum member mass
is not strikingly different from the size-of-sample hypothesis (Larson, 1982),
so it is not trivially dismissed. Distinct high-mass truncations are evident
in the raw present-day mass functions of 30 Dor and Carina OB2 (Weidner &
Kroupa, 2004; Wright et al., 2010), but with the severe caveat that the
present-day mass function of these $>4\,\rm Myr$ old systems will also be
truncated from the death of the earliest-forming massive stars. This issue of
sensitivity to the uncertain star formation history (among others) muddles the
interpretation of observations.

The situation in the simulations is clearer. In §3.1 we showed that our 100
low-mass cluster formation simulations, taken together with previous STARFORGE
simulations of other conditions, suggest an environmentally-varying high-mass
truncation in the IMF that is generally distinct from the usually-assumed
universal $100-150M_{\odot}$ cutoff. We do not yet have sufficient coverage of
parameter space to map out the detailed scaling, but the $10\times$ more-
massive cloud IMF shown in Figure 1 indicates that at least total cloud mass
is a significant factor.

One factor affecting the maximum stellar mass in a given cloud could simply be
the time that a star has to draw gas from the cloud. These simulated clouds
tend to disperse due to feedback after a time proportional to the free-fall
timescale $t_{\rm ff}$, and at fixed mean surface density $t_{\rm ff}\propto
M^{1/4}$, so more-massive clouds tend to live longer, and more-extended
accretion is possible. At first glance this explanation seems tempting for the
cases shown in Fig. 1, because the $\sim 10\times$ more-massive cloud produces
a star $\sim 10^{1/4}=1.8$ times more massive.

But the cloud lifetime alone does not explain all variation, and predicts the
wrong trend in other cases: Guszejnov et al. (2022) found their densest,
shortest-lived cloud produced more-massive stars than the fiducial case. The
other facet of the problem is feedback: simulations of these same clouds
without feedback were easily able to produce 100+$M_{\odot}$ stars (Guszejnov
et al., 2020; Guszejnov et al., 2021). The relation between the cutoff and
extensive cloud properties may therefore be understood in terms of the way
massive stars tend to form in these simulations: from accretion of extended
gas flows over timescales as long as several Myr (Grudić et al., 2022), which
are hardly separated in length- or time-scale from the parent cloud size and
crossing time. The accretion of massive stars and the resulting influence of
feedback on cloud dynamics are therefore coupled. Denser clouds are generally
more resistant to feedback, and may therefore produce more-massive stars if
all else is equal.

####  4.1.1 Analytic model

To build some analytic intuition for how feedback and cloud properties might
influence the upper cutoff of the IMF, let us consider how the stellar mass
required to disrupt a cloud might scale under some simplifying assumptions. We
estimate an upper bound on $M_{\rm up}$ by determining the mass that a star
must have to single-handedly disrupt a cloud within some fraction of a
freefall time $N_{\rm ff}\lesssim 1$, i.e. soon enough to prevent significant
additional star formation.

If the cloud is disrupted by ionizing radiation, we may adopt the Hosokawa &
Inutsuka (2006) D-type HII region expansion solution and ask what H ionizing
photon production rate $\mathcal{Q}\left(H^{\rm 0}\right)$ is required for the
bubble to sweep out the cloud radius $R$ within $N_{\rm ff}$ freefall times,
$t_{\rm ff}=\sqrt{3\uppi/\left(32G\rho\right)}$. Assuming $R$ is much larger
than the Strömgren radius, we obtain

$\mathcal{Q}\left(H^{\rm 0}\right)\gtrsim\frac{6912G^{2}M^{3/2}m_{\rm
p}^{2}X_{\rm H}^{2}\alpha_{\rm B}\left(T_{\rm i}\right)\mu_{\rm
i}^{2}\Sigma^{5/2}}{2401\uppi^{5/2}N_{\rm ff}^{4}k_{\rm B}^{2}T_{\rm i}^{2}},$
(2)

where $M$ is the mass of the cloud, $X_{\rm H}$ the mass fraction of H,
$\mu_{\rm i}\sim 0.61$ the mean molecular weight of ionized gas,
$\Sigma=M/\uppi R^{2}$ the mean cloud surface density, $T_{\rm i}\sim
10^{4}\rm K$ the temperature of ionized gas, and $\alpha_{\rm B}\left(T_{\rm
i}\right)\approx 2.6\times 10^{-13}\left(T_{\rm i}/10^{4}\rm
K\right)^{-0.85}\,\rm cm^{3}\,s^{-1}$ is the case B recombination coefficient
(Osterbrock & Ferland, 2006, Table 2.1). More conveniently,

$\mathcal{Q}\left(H^{\rm 0}\right)\gtrsim 10^{47}\,\mathrm{s}^{-1}N_{\rm
ff}^{-4}\Sigma_{\rm 2}^{5/2}M_{\rm 4}^{3/2}T_{\rm i,4}^{-2.85},$ (3)

where $\Sigma_{\rm 2}=\Sigma/10^{2}M_{\odot}\rm pc^{-2}$, $M_{\rm
4}=M/10^{4}M_{\odot}$, and $T_{\rm i,4}=T/10^{4}\rm K$. Approximating the
Diaz-Miller et al. (1998) scaling of ionizing flux with stellar mass222In
identifying $\mathcal{Q}\left(H^{\rm 0}\right)$ with a certain stellar mass we
have assumed that all of the ionizing flux originates in the most massive
star. In practice, the most massive star accounts for $>50\%$ of the total
ionizing flux in 95 of our 100 simulations, so this is not an especially bad
assumption. In the stellar mass range $\sim 10-30M_{\odot}$
$\mathcal{Q}\left(H^{\rm 0}\right)$ scales so steeply with stellar mass that
the most massive star tends to dominate. over the range $10-30M_{\odot}$ as a
power-law $\log\mathcal{Q}\left(H^{\rm
0}\right)=44.7+7.3\log\left(M_{\star}/10M_{\odot}\right)$, we find

$M_{\star}\gtrsim 16M_{\odot}N_{\rm ff}^{-0.55}\Sigma_{\rm 2}^{0.34}M_{\rm
4}^{0.2}T_{\rm i,4}^{-0.39}$ (4)

Hence the maximum stellar mass scales within increasing surface density and
cloud mass. For the 100 simulations with $M_{\rm 4}=0.2$ and $\Sigma_{\rm
2}=0.7$, the typical maximum cluster mass was $\sim 18M_{\odot}$, implying
$N_{\rm ff}\sim 0.4$, fairly consistent with the cloud disruption timescales
documented for these simulations (Grudić et al., 2022; Guszejnov et al.,
2022). Also notable is the dual role of $T_{\rm i}$, setting both the gas
pressure driving bubble expansion and the efficiency of recombination.
Dramatically lower ionizing fluxes are required to disrupt a cloud if $T_{\rm
i}$ is greater, e.g. in very low-metallicity clouds lacking metal line
coolants. This agrees qualitatively with the $0.01Z_{\odot}$ simulations
presented in Guszejnov et al. (2022), which did generally have significantly
lower maximum stellar masses (and star formation efficiency) than their
$Z_{\odot}$ counterparts.

Photoionization feedback is not universally effective across all of cloud
parameter space, so other mechanisms could also affect the maximum stellar
mass. The picture may change qualitatively when scaling to sufficiently
massive (e.g. $\gtrsim 10^{6}M_{\odot}$) clouds, with freefall times exceeding
$10\rm Myr$, both because the greater escape speed relative to the ionized gas
sound speed will make a given amount of feedback less efficient, and because
massive stars can die over the longer timescale.

It is possible to repeat this exercise substituting feedback bubble solutions
appropriate for other feedback mechanisms single-scattered radiation pressure
(Steigman et al., 1975; Krumholz & Matzner, 2009; Raskutti et al., 2016) and
stellar winds (Weaver et al., 1977; Lancaster et al., 2021). Those feedback
mechanisms also scale very steeply with stellar mass, so qualitatively-similar
scaling laws for the required stellar mass will generally result. Those
feedback mechanisms do not contribute significantly to cloud disruption in the
simulations discussed here, so no meaningful comparison to model predictions
can be made. As stars approach the Eddington limit ($\gtrsim 30M_{\odot}$)
multiple-scattered infrared radiation pressure should also become important
(Larson & Starrfield, 1971; Wolfire & Cassinelli, 1987; Kuiper et al., 2010),
likely tapering off the scaling law in a manner that is sensitive to dust
properties. But how the bound on the maximum stellar mass from IR radiation
pressure connects to the large-scale cloud properties is less clear, because
it operates mainly on much smaller scales: within the dense protostellar
envelope, just beyond the dust destruction front at a radius of a few AU.

####  4.1.2 Broader implications

If different cloud properties imprint different upper truncations in the IMF,
then there are important implications for the stellar populations of galaxies.
GMCs span a broad range of masses and densities within even a single galaxy
(Heyer & Dame, 2015; Freeman et al., 2017; Sun et al., 2018), so the GMC-level
IMFs would vary in their upper truncation and the resulting galaxy-wide IMF
would have a different shape. In particular, if $M_{\rm up}$ correlates
positively with the host cloud mass then the galactic IMF will generally be
steeper than the individual cluster IMFs. This is one way in which our models
may be able to reproduce the observed high-mass slope in the field population,
which is probably steeper than our measured $\alpha=-2$. This is also a key
prediction of the integrated galactic IMF (IGIMF) theory developed by Weidner
& Kroupa (2006) and many subsequent works (see Kroupa et al. 2013 for review).
However, our eventual star cluster demographics are more random than the
“optimal sampling” model usually assumed in such works. Crucially, the
simulations have a scatter of $0.13\,\rm dex$ in the relation between cluster
mass and the most-massive star, whereas optimal sampling has none. This
variation compounds to the $\sim 0.4\,\rm dex$ of scatter in luminosity shown
in Figure 3. Nevertheless, it is the systematic variation in $M_{\rm up}$ –
and not so much the specifics of the sampling process – that is required for
the composite IMF to have a different shape from the cluster-level IMF, and
our simulations give physical motivation to this hypothesis. Whether the
corrections to normal population synthesis are large or small in practice
remains to be seen, as it will depend on how specifically $M_{\rm up}$ varies
across the vast parameter space of star-forming cloud properties.

Because the fully-formed clusters’ photometric properties are similar to
stochastic sampling from a truncated IMF in the simulations, our result is
similar to the one tested observationally by Andrews et al. (2013). They
compared photometric models stochastically sampled from a Kroupa (2001) IMF
assuming $M_{\rm up}=30M_{\odot}$ and $120M_{\odot}$ respectively, finding
that $M_{\rm up}=120M_{\odot}$ is better able to account for the photometry of
$\sim 10^{3}M_{\odot}$ young star clusters in NGC 4214. However, we emphasize
that 1. our clusters are an order of magnitude less massive, and 2. the
simulations suggest that the $M_{\rm up}$ should vary with the properties of
the cluster or its progenitor cloud, and $\alpha$ can also vary from the exact
Salpeter (1955) value. Thus it is not surprising that a model with a hard
cutoff of $30M_{\odot}$ with $\alpha=-2.35$ exactly across the entire galaxy
is in tension with the data. Given the widely-varying properties of giant
molecular clouds, we expect there should be conditions in most galaxies
capable of producing stars more massive than $30M_{\odot}$.

In principle, star cluster photometry may constrain the variation of $\alpha$
and/or $M_{\rm up}$ in greater detail, e.g. by incorporating a model for their
scaling and intrinsic scatter in a full Bayesian forward-model (e.g. Krumholz
et al., 2019). However in practice $\alpha$ and $M_{\rm up}$ tend to be
degenerate because both parameters mainly affect the ratio of ionizing (and
hence $\rm H\alpha$) luminosity to bolometric luminosity. Photometric models
allowing even $\alpha$ alone to vary freely tend to be degenerate (Ashworth et
al., 2017), so simulations may be useful for establishing physically-motivated
priors to break degeneracies.

###  4.2 Implications for simulations and feedback modeling

Stochastic sampling from a fixed standard IMF has become increasingly popular
in numerical simulations of galaxies and star cluster formation, with the aim
of modeling feedback, stellar dynamics, and chemical enrichment more
accurately on small scales. But in §3.1 we found that the conditions we
simulated do not allow the formation of $\gtrsim 30M_{\odot}$ stars, requiring
an environmentally-dependent upper IMF truncation to account for the range of
stellar masses observed. Ours is likely an extreme case for very low-mass
clusters, but nevertheless the variation of $M_{\rm up}$ could be important
for feedback modeling in general. Fortunately the specific choice of $M_{\rm
up}$ will not have a large effect on the rate of core-collapse supernovae from
a given stellar population, because the number of progenitors is dominated
by$\sim 8M_{\odot}$ stars, well below the cutoff seen even in these low-mass
clusters. Therefore the overall energy budget of supernova feedback will not
be highly sensitive.

But a lower value of $M_{\rm up}$ would produce significantly less radiative
feedback from photoionization and early-time radiation pressure. Compared to a
Kroupa (2001) IMF with $\alpha=-2.3$ and $M_{\rm up}=150M_{\odot}$, a cluster
with our sample IMF with $\alpha=-2$ and $M_{\rm up}\sim 28M_{\odot}$ would
have a zero-age bolometric luminosity $5\times$ less, and produce $10\times$
fewer H-ionizing photons. There would also be significantly fewer stars
producing long-lived, powerful Wolf-Rayet winds (Meynet & Maeder, 2005). Since
these early feedback mechanisms are key for regulating star formation on cloud
and star cluster scales, cloud and star cluster properties and the pre-
processing of the ISM for supernova explosions may be affected. A deficit of
very massive stars would also result in measurably-different chemical
abundances from supernova ejecta (e.g. McWilliam et al., 2013).

Even if the correct IMF is known, we also found in §3.2 that the way in which
stars of different masses are sampled as the cluster formed is not well-
modeled by random sampling. The discrepancy between the number of currently-
massive stars and random sampling is partly due to the non-negligible time
required for massive stars to form. Hence it is possible that stochastic
models could be augmented with a model accounting for the finite time required
to assemble a star. Without accounting for the time delay, massive stars would
be predicted to form too early in a cluster’s lifetime, and their feedback
would influence subsequent star formation disproportionately.

###  4.3 Caveats

Star formation simulations require many physical approximations and can only
treat a segment of the physical dynamic range self-consistently, so their
results should be interpreted carefully. A detailed discussion of known and
potential caveats and uncertainties of the STARFORGE numerical model is given
in Grudić et al. (2021) §5.2.

A potentially-important uncertainty is the initial and boundary conditions.
Here we simulate an ensemble of isolated, initially-uniform, spherical clouds
with Gaussian random initial turbulence, and real GMCs are unlikely to be
well-described by any of these properties. Some studies have found that such
spherical models may exhibit realistic properties and star formation dynamics
once clumpy substructure has had time to form (Rey-Raposo et al., 2015;
Priestley et al., 2023). But any dense cloud must inevitably have some
continuity with the converging galactic gas flow that formed it, so isolated
cloud evolution is unlikely to capture the full dynamics of star formation.
The problem of the IMF and its sampling process should be revisited with
alternative initial/boundary conditions setups (e.g., Lane et al., 2022), and
ideally a larger-scale setup embedded within a galaxy.

And lastly, we re-emphasize that we have intentionally simulated a very
particular case of low-mass star formation to stress-test the usual
assumptions of stellar population synthesis, and that this is not
representative of the massive complexes that form most stars. We have shown
that discrepancies with stochastic sampling exist in detail, and this should
not change qualitatively with the scale of the star-forming region. But it
remains to be seen whether the differences from stochastic sampling are large
in general. For instance, if the sampling always regresses to the mean once
$\sim 100$ stars have formed then the correction would be negligible in any
massive star-forming complex. But 100 just happens to be the typical cluster
size here, so it could also be that clusters regress to the mean on Fig. 2
only once most stars have formed, in which case the difference from stochastic
sampling will be important at all scales. It will be possible to clarify this
issue with an upcoming statistics suite of larger clouds, currently in
progress.

##  5 Conclusion

In this work we use numerical simulations to test two common assumptions when
modeling stellar populations: that stellar masses are randomly, independently
sampled, and that the underling IMF is invariant. Within our models, both
assumptions break down to some extent:

  * •

There is a high-mass truncation in the IMF $M_{\rm up}$ that can generally be
significantly lower than largest known stellar masses (Fig. 1), at least in
low-mass clusters. This truncation varies with cloud properties such as mass,
so different environmental conditions are required to account for the full
observed range of stellar masses, and in general it is not valid to assume
that all clouds can form arbitrarily-massive stars.

  * •

Stars of different masses do not form in a random order. Massive stars in the
simulated clouds begin forming earlier and finish forming later than average,
so assuming random sampling during star formation overestimates the number of
massive stars that exist at a given time.

Although star formation is not random and uncorrelated in our deterministic
simulations, stochastic sampling does appear to provide an adequate
description of the end result of star formation, with the major caveat that
the high-mass slope and truncation for the system in question must be known,
or generally allowed to vary as a model parameter. Stars do not form through a
roll of the dice, but the laws of physics may conspire to make it appear so.

## Acknowledgements

We thank Zhiqiang Yan, Bruce Elmegreen, Martijn Wilhelm, Daniela Calzetti, and
Pavel Kroupa for enlightening exchanges that motivated this work. Support for
MYG was provided by NASA through the NASA Hubble Fellowship grant #HST-
HF2-51479 awarded by the Space Telescope Science Institute, which is operated
by the Association of Universities for Research in Astronomy, Inc., for NASA,
under contract NAS5-26555. SSRO was supported by NSF through CAREER award
1748571, AST-2107340 and AST-2107942, by NASA through grants 80NSSC20K0507 and
80NSSC23K0476 and by the Oden Institute through a Moncrief Grand Challenge
award. CAFG was supported by NSF through grants AST-2108230 and CAREER award
AST-1652522; by NASA through grants 17-ATP17-0067 and 21-ATP21-0036; by STScI
through grant HST-GO-16730.016-A; by CXO through grant TM2-23005X; and by the
Research Corporation for Science Advancement through a Cottrell Scholar Award.
Support for DG was provided by NASA through the NASA Hubble Fellowship grant
#HST-HF2-51506 awarded by the Space Telescope Science Institute. Support for
PFH was provided by NSF Research Grants 1911233, 20009234, 2108318, NSF CAREER
grant 1455342, NASA grants 80NSSC18K0562, HST-AR-15800. The suite of 100 new
simulations was run on the Anvil supercomputer, through XSEDE award TG-
PHY210070. Analysis was performed using Frontera allocation AST21002. This
research is part of the Frontera computing project at the Texas Advanced
Computing Center. Frontera is made possible by National Science Foundation
award OAC-1818253.

## Data Availability

The complete dataset of simulation snapshots with time-dependendent stellar
properties is available at this link. A public version of the GIZMO code is
available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.

## References

  * Andrews et al. (2013) Andrews J. E., et al., 2013, ApJ, 767, 51
  * Ashworth et al. (2017) Ashworth G., et al., 2017, MNRAS, 469, 2464
  * Bastian et al. (2010) Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339
  * Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Chabrier (2005) Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library Vol. 327, The Initial Mass Function 50 Years Later. p. 41
  * Crowther et al. (2016) Crowther P. A., et al., 2016, MNRAS, 458, 624
  * Diaz-Miller et al. (1998) Diaz-Miller R. I., Franco J., Shore S. N., 1998, ApJ, 501, 192
  * Dinnbier & Walch (2020) Dinnbier F., Walch S., 2020, MNRAS, 499, 748
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Freeman et al. (2017) Freeman P., Rosolowsky E., Kruijssen J. M. D., Bastian N., Adamo A., 2017, MNRAS, 468, 1769
  * Fujii & Portegies Zwart (2015) Fujii M. S., Portegies Zwart S., 2015, MNRAS, 449, 726
  * Fumagalli et al. (2011) Fumagalli M., da Silva R. L., Krumholz M. R., 2011, ApJ, 741, L26
  * Grudić & Hopkins (2019) Grudić M. Y., Hopkins P. F., 2019, MNRAS, 488, 2970
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216
  * Guszejnov et al. (2020) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2020, MNRAS, 496, 5072
  * Guszejnov et al. (2021) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 502, 3646
  * Guszejnov et al. (2022) Guszejnov D., Grudić M. Y., Offner S. S. R., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2022, MNRAS, 515, 4929
  * Heyer & Dame (2015) Heyer M., Dame T., 2015, Annual Review of Astronomy and Astrophysics, 53, 583
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins (2016) Hopkins P. F., 2016, MNRAS, 462, 576
  * Hopkins (2018) Hopkins A. M., 2018, PASA, 35, e039
  * Hopkins & Grudić (2019) Hopkins P. F., Grudić M. Y., 2019, MNRAS, 483, 4187
  * Hopkins & Raives (2016) Hopkins P. F., Raives M. J., 2016, MNRAS, 455, 51
  * Hopkins et al. (2020) Hopkins P. F., Grudić M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MNRAS, 491, 3702
  * Hopkins et al. (2023) Hopkins P. F., et al., 2023, MNRAS, 519, 3154
  * Hosokawa & Inutsuka (2006) Hosokawa T., Inutsuka S.-i., 2006, ApJ, 648, L131
  * Hu et al. (2017) Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa (2002) Kroupa P., 2002, Science, 295, 82
  * Kroupa & Jerabkova (2021) Kroupa P., Jerabkova T., 2021, arXiv e-prints, p. arXiv:2112.10788
  * Kroupa et al. (2013) Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, The Stellar and Sub-Stellar Initial Mass Function of Simple and Composite Populations. p. 115, doi:10.1007/978-94-007-5612-0_4
  * Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, The Astrophysical Journal, 703, 1352
  * Krumholz et al. (2015) Krumholz M. R., et al., 2015, ApJ, 812, 147
  * Krumholz et al. (2019) Krumholz M. R., Adamo A., Fumagalli M., Calzetti D., 2019, MNRAS, 482, 3550
  * Kuiper et al. (2010) Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ, 722, 1556
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2
  * Lancaster et al. (2021) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021, ApJ, 914, 90
  * Lane et al. (2022) Lane H. B., Grudić M. Y., Guszejnov D., Offner S. S. R., Faucher-Giguère C.-A., Rosen A. L., 2022, MNRAS, 510, 4767
  * Larson (1982) Larson R. B., 1982, MNRAS, 200, 159
  * Larson & Starrfield (1971) Larson R. B., Starrfield S., 1971, A&A, 13, 190
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Lewis et al. (2023) Lewis S. C., et al., 2023, ApJ, 944, 211
  * Liow et al. (2022) Liow K. Y., Rieder S., Dobbs C. L., Jaffa S. E., 2022, MNRAS, 510, 2657
  * McWilliam et al. (2013) McWilliam A., Wallerstein G., Mottini M., 2013, ApJ, 778, 149
  * Meynet & Maeder (2005) Meynet G., Maeder A., 2005, A&A, 429, 581
  * Offner et al. (2009) Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ, 703, 131
  * Offner et al. (2014) Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, Protostars and Planets VI, pp 53–75
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Priestley et al. (2023) Priestley F. D., Clark P. C., Whitworth A. P., 2023, MNRAS, 519, 6392
  * Raskutti et al. (2016) Raskutti S., Ostriker E. C., Skinner M. A., 2016, ApJ, 829, 130
  * Rey-Raposo et al. (2015) Rey-Raposo R., Dobbs C., Duarte-Cabral A., 2015, MNRAS, 446, L46
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Schechter (1976) Schechter P., 1976, ApJ, 203, 297
  * Schwarz (1978) Schwarz G., 1978, Annals of Statistics, 6, 461
  * Smith (2021) Smith M. C., 2021, MNRAS, 502, 5417
  * Steigman et al. (1975) Steigman G., Strittmatter P. A., Williams R. E., 1975, ApJ, 198, 575
  * Su et al. (2018) Su K.-Y., et al., 2018, MNRAS, 480, 1666
  * Sun et al. (2018) Sun J., et al., 2018, ApJ, 860, 172
  * Tang et al. (2023) Tang J., Grasha K., Krumholz M. R., 2023, arXiv e-prints, p. arXiv:2301.05912
  * Tout et al. (1996) Tout C. A., Pols O. R., Eggleton P. P., Han Z., 1996, MNRAS, 281, 257
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Weidner & Kroupa (2004) Weidner C., Kroupa P., 2004, MNRAS, 348, 187
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Wolfire & Cassinelli (1987) Wolfire M. G., Cassinelli J. P., 1987, ApJ, 319, 850
  * Wright et al. (2010) Wright N. J., Drake J. J., Drew J. E., Vink J. S., 2010, ApJ, 713, 871
  * da Silva et al. (2012) da Silva R. L., Fumagalli M., Krumholz M., 2012, ApJ, 745, 145