# Cluster Population Demographics in NGC 628 Derived from Stochastic
Population Synthesis Models

Jianling Tang,1 Kathryn Grasha,1,2 Mark R. Krumholz1,2
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The physical properties of star cluster populations offer valuable insights
into their birth, evolution, and disruption. However, individual stars in
clusters beyond the nearest neighbours of the Milky Way are unresolved,
forcing analyses of star cluster demographics to rely on integrated light, a
process that is fraught with uncertainty. Here infer the demographics of the
cluster population in the benchmark galaxy NGC 628 using data from the Legacy
Extra-galactic UV Survey (LEGUS) coupled to a novel Bayesian forward-modelling
technique. Our method allows analysis of a total of 1178 clusters in the LEGUS
catalogue, roughly a factor of 4 more than previous studies that required
severe completeness cuts to the data. Our results indicate that the cluster
mass function is truncated at $\sim 10^{4}$ $\mathrm{M}_{\odot}$, consistent
with proposed relations between truncation mass and star formation surface
density. We find that cluster disruption begins early, $\sim 10$ Myr after
formation, but that it is relatively mild, with clusters requiring on average
$2-3$ times their present age to disrupt; we do not find any evidence for mass
dependent disruption. We also do not find convincing evidence for any radial
variations in these conclusions, though we find suggestive hints that inner
galaxy clusters may be more prone to disruption. Confirming or refuting these
hints will require future observations to increase the sample size of outer
galaxy clusters.

###### keywords:

galaxies: clusters: general – techniques: photometric – software: data
analysis – methods: statistical

††pubyear: 2023††pagerange: Cluster Population Demographics in NGC 628 Derived
from Stochastic Population Synthesis Models–A

##  1 Introduction

The distribution of star clusters’ properties (mass, age, physical size) and
the variation of this distribution with galactic environment provides crucial
clues to the physics of star formation and star cluster evolution. Moreover,
because of the processes responsible for setting this distribution depend on
galactic environment, cluster demographics can also trace the history of
galaxy assembly and evolution, acting as a fossil record of the environments
that existed when clusters formed. Given the importance of cluster
demographics, it is not surprising that there have been many attempts to
measure them, and that cluster demographics figure prominently in the science
cases for a number of large surveys of nearby galaxies, such as PHAT
(Dalcanton et al., 2012), PHANGS-HST (Lee et al., 2022), and LEGUS (Calzetti
et al., 2015a); see Krumholz et al. (2019a) for a detailed review.

Different aspects of cluster demographics probe different physics. By studying
the shape of the cluster mass function (CMF), and in particular the initial
CMF (ICMF) that applies to the youngest clusters, we can constrain star
formation theories and provide an observational check on simulations. Studies
to date show that over much of its range the ICMF is well-described by a
powerlaw $dN/dM\propto M^{\alpha_{M}}$ with a slope $\alpha_{M}\approx-2$
across a wide range of galaxy properties, corresponding to equal mass per
logarithmic bin and suggesting a scale-free formation process (e.g., Fall &
Chandar, 2012; Chandar et al., 2014). However, the shape of the ICMF at its
high-mass end ($\gtrapprox 10^{4}M_{\odot}$) remains uncertain. Some authors
report that a truncated distribution such as a Schechter function describes
the data better than the pure power-law (e.g., Bastian et al., 2012; Adamo et
al., 2017; Johnson et al., 2017), while others question the statistical
significance of claimed detections and instead suggest that the dearth of
massive clusters is simply a size-of-sample effect (e.g., Larsen, 2002; Mok et
al., 2020). If there is a real truncation in the mass function, its location
must depend somehow on the galactic environment, since rapidly star-forming
galaxies with large cluster populations often harbour clusters with masses
larger than the reported truncation masses in more modestly star-forming
galaxies (e.g., Chandar et al., 2015; Linden et al., 2017). The presence of a
non-existence of a high-mass break in the ICMF, its value, and its variation
with galactic environment provide an important clue to the star formation
process.

As clusters evolve and disperse from their birthplaces, they experience mass
loss. Various theories govern how they disrupt, ranging from the “infant
mortality” stage of rapid gas removal after star formation (Hills, 1980;
Baumgardt & Kroupa, 2007; Dinnbier & Walch, 2020) to long-timescale processes
such as two-body relaxation (Fall & Zhang, 2001; Lamers et al., 2010a;
Kruijssen et al., 2012) to processes such as external tidal shocking that
operate on intermediate timescales (Lamers et al., 2005b; Bastian et al.,
2005; Gieles et al., 2006; Gieles & Renaud, 2016; Howard et al., 2017; Fensch
et al., 2019; Webb et al., 2019). In order to place observational constraints
on these theories, we study the cluster age function (CAF). As with the ICMF,
observations suggest that the CAF is reasonably well-described by a powerlaw
form $dN/dT\propto T^{\alpha_{T}}$, where $T$ is cluster age. However, the
value of $\alpha_{T}$ remains highly controversial, at least in part due to
disagreement among observational groups about what constitutes a cluster
(Krumholz et al., 2019a). Despite this uncertainty, the CAF encodes the
timescales of diverse cluster disruptive processes. One question of particular
interest is the form of the joint mass-age distribution, which gives a
complete description of cluster formation and disruption. If the mass and age
distributions are separable, i.e. if $d^{2}N/dM\,dT\propto(dN/dM)(dN/dT)$,
this implies mass-independent disruption (MID) of clusters (Fall et al., 2005,
2009; Chandar et al., 2015), while if they are not this implies mass-dependent
disruption (MDD; Lamers et al. 2005a; Gieles et al. 2007; Grudić et al. 2021).
Knowing whether MID or MDD holds in a particular galaxy or sub-galactic region
would in turn place strong constraints on the mechanisms by which clusters
disrupt.

Part of the reason that observations have not yet resolved questions such as
the existence and value of a cutoff in the ICMF, the slope of the CAF, and
whether cluster disruption is mass-dependent or -independent, is that
determining the properties of clusters from observations is not trivial. In
order to sample a wide range of environments, studies of cluster populations
must work with integrated light rather than resolved stellar populations,
since in the crowded environments of star clusters it is generally only
possible to resolve individual stellar sources for a handful of the most
nearby galaxies. The traditional approach to extracting cluster demographics
from this type of data is to convert the integrated light measurements for
each cluster to physical properties such as mass and age by comparing the
observed photometry to a grid of simple stellar population (SSP) models and
generating a set of best-fitting parameters. One then bins the clusters by
mass and age to obtain mass and age distributions.

However, this approach encounters several difficulties. First, the binning
process usually entails the loss of useful information, and as a result
parameters determined by fitting data to histograms are heavily biased
(Maschberger & Kroupa, 2009; El-Badry et al., 2017). Second, the conventional
approach of using $\chi^{2}$ fitting to convert photometry to masses and ages
implicitly assumes that the uncertainties on cluster mass and age can be
approximated as Gaussian. This is often a poor assumption, because the mapping
between physical properties and photometry is non-linear and non-monotonic,
particularly once one adds the additional complication of dust extinction.
Consequently, a particular set of photometric measurements may plausibly fit
two (or more) widely-separated loci in physical parameter space, yielding
posterior mass and age distributions with complex non-symmetrical and multi-
modal shapes (de Meulenaer et al., 2014; Fouesneau et al., 2014; Krumholz et
al., 2015b; Bialopetravičius et al., 2019). The problem is exacerbated for
low-mass clusters, where stochastic sampling of the IMF leads to a wide range
of possible photometry even for clusters of fixed mass and age (Piskunov et
al., 2009).

However, the most severe issue for measuring cluster demographics is
completeness. The problem is that star clusters become fainter as one moves to
both lower mass and older age, so a magnitude limit corresponds to a complex
shape in the parameter space of age and mass, meaning that both the mass and
age distributions are subject to large potential biases. The conventional way
of handling this is to discard clusters outside of a limited range of both
mass and age over which completeness is expected to be $\gtrsim 90\%$, e.g.,
to retain only clusters with estimated masses $5,000$ M⊙ and estimated ages
$<200$ Myr. Such a drastic truncation of the sample avoids bias, but at the
cost of a large loss of statistical power at both lower masses and old ages.
The former is particularly concerning, since the steep mass function means
that low-mass clusters form the majority of the available sample; the sample
truncation required to avoid bias from survey magnitude limits therefore often
involves discarding the majority of the data.

The main objective of this paper is to present a complete analysis pipeline
and preliminary results for cluster demographics in the benchmark galaxy NGC
628 using the catalogue of clusters measured by the _Hubble Space Telescope_
Treasury Program Legacy Extragalactic Ultraviolet Survey (LEGUS; Calzetti et
al., 2015a) coupled to the novel Bayesian forward modelling method proposed by
Krumholz et al. (2019b) that overcomes the problems identified above.
Specifically, this method naturally copes with complex and non-deterministic
mappings between photometric measurements and physical properties, and it
enables us to use a total of 1178 clusters catalogued by LEGUS with non-zero
completeness values, as compared to earlier modelling where severe
completeness cuts reduced the sample to $\approx 300$ (Adamo et al., 2017).
This paper is structured as follows: in Section 2, we introduce our target
galaxy, cluster catalogue, and analysis of completeness. In Section 3, we
summarise our analysis method. We present our results for cluster demographics
in Section 4, and discuss their implications in Section 5. We summarise the
findings and discuss future prospects in Section 4.

##  2 Data Description

Here we summarise the properties of our target galaxy (Section 2.1), the star
cluster catalogue we use as the basis for our study (Section 2.2), and our
treatment of survey completeness (Section 2.3).

###  2.1 NGC 628

Our target NGC 628, also known as Messier 74, is a well-studied grand-design
spiral galaxy, located at a distance of 9.9 Mpc (Olivares E. et al., 2010). We
choose NGC 628 for this study as its large size (apparent radius of 5.2’;
Gusev et al. 2014) and face-on orientation allow for a detailed study of the
stellar cluster population. We use observations of NGC 268 taken from the
LEGUS survey, a Cycle 21 HST Treasury programme that targeted 50 local
galaxies ($\lesssim$15 Mpc) with the Hubble Space Telescope with broad band
filter coverage from the UV to the near IR. All targets were imaged with
either the Wide Field Camera 3 (WFC3) or the Advanced Camera for Surveys (ACS)
in the NUV (WFC3 F275W), U-band (F336W), B-band (ACS/F435W or WFC3/F438W),
V-band (ACS or WFC3/F555W or F606W), and I-band (ACS or WFC3 F814W). For the
remainder of this paper we follow the conventional Johnson passband naming
convention UV, U, B, V and I without converting to the Johnson system. As part
of the LEGUS survey, NGC 628 is covered by two pointings; one at the galactic
centre (NGC 628c) and the other at the eastern edge (NGC 628e). We show the
LEGUS V-band image of NGC 628 in Figure 1. Detailed descriptions of the survey
and data reduction are provided in Calzetti et al. (2015b).

Figure 1: V-band (F555W) mosaic image of the centre and east fields of NGC 628
from LEGUS-HST. Symbols show the location of the star clusters, colored by the
cluster classifications (Section 2.2) : class 1 (red circle), class 2 (blue
circle), and class 3 (blue circle). For detailed description of cluster
classification method, Section 2.2.

###  2.2 Cluster Catalogue

The photometric star cluster catalogue for NGC 628 that we use in this study
comes from the LEGUS survey, and we refer readers to Grasha et al. (2015) and
Adamo et al. (2017) for a complete description of how it is constructed. Here,
we provide a brief description of the parts of the process most pertinent to
this study. LEGUS uses a two-step pipeline for cluster classification. First,
an automated cluster extraction tool is used to extract potential cluster
candidates with at least a 3$\sigma$ detection in a minimum of 5 contiguous
pixels from the white light images. The cluster candidates are then further
refined and required to have a V-band concentration index (CI; the difference
in magnitude between a radius of 1 and 3 pixels) of $>1.4$ mag for the centre
pointing and $>1.3$ mag for the east pointing. The CI cut serves to separate
point-like sources (i.e., stars) from more extended sources (i.e., star
clusters). The star cluster candidates are required to be detected in the
V-band and a minimum of at least three of the other filters, with photometric
error $\sigma_{\lambda}\leq 0.3$ mag in each band.

To create the final cluster catalog, the LEGUS team then visually inspects all
candidates brighter than $-$6 in the V-band and classifies them into four
morphological categories: (1) centrally concentrated clusters with spherically
symmetric profiles; (2) clusters with asymmetric radial profiles; (3) multi-
peaked clusters with underlying diffuse emission; and (4) non-cluster
contaminants, such as background galaxies, stars, bad pixels, or edge
artefacts. In this work, we limit our analysis to Class 1 (426), 2 (437), and
3 (413) clusters, for a total sample of 1276; we reduce this to 1275 by
removing one cluster that we found to be a duplicate entry in the LEGUS NGC
628c and 628e catalogues, located in the small region where the fields
overlap. Figure 1 shows the mosaic V-band image of NGC 628 overlaid with the
cluster positions, coloured by their morphological classification.

###  2.3 Observational Completeness Limits

Because we intend to forward model the cluster population, we require
knowledge of the completeness of the observational catalogue. Specifically, we
require knowledge of the function $P_{\mathrm{obs}}(\mathbf{m})$, which
describes the probability that a hypothetical star cluster with a vector of
magnitudes $\mathbf{m}$ in the various LEGUS filters would be included in the
cluster catalogue. To compute this function, we employ the completeness limits
reported by Adamo et al. (2017), who carry out artificial cluster tests to
derive the completeness of the LEGUS automated catalogue generation procedure
in each filter independently. For each filter $F$ and a range of magnitudes
$m_{F}$ in that filter, they determine the probability
$P_{\mathrm{obs},F}(m_{F})$ that a cluster would be recovered by the automated
extraction procedure. In what follows we linearly interpolate these tabulated
data to obtain a continuous function $P_{\mathrm{obs},F}(m_{F})$ that gives
the probability that a cluster of arbitrary magnitude $m_{F}$ in filter $F$
will be recovered.

To determine the completeness for a cluster with a vector of magnitudes
$\mathbf{m}$, we use 10,000 Monte Carlo trials. In each trial, we randomly
assign each filter $F$ to be a detection or a non-detection with probability
$P_{\mathrm{obs},F}(m_{F})$ as determined from the interpolated artificial
cluster test results. We then determine from this set of detections and non-
detections if the cluster would be catalogued following the same criteria used
in construction of the actual LEGUS catalogue, i.e., the cluster is catalogued
only if it is (1) detected in V-band and at least three other bands, and (2)
has a visual magnitude $M_{\mathrm{V}}\leq-6$. We then take
$P_{\mathrm{obs}}(\mathbf{m})$ to be equal to the fraction of the Monte Carlo
trials in which the cluster is catalogued.

We use this method to calculate the completeness both of our library of
synthetic clusters (see Section 3) and of the actual LEGUS catalogue. For the
latter, we find a total of 97 clusters for which we estimate
$P_{\mathrm{obs},F}(\mathbf{m})=0$; these can be present in the catalogue
because a small number of clusters with $m_{\mathrm{V}}>-6$ were added by
hand. To avoid introducing errors in our completeness estimate we remove these
from the sample. We remove a total of 97 clusters from the sample of NGC 628
due to their zero completeness values, resulting in a final sample size of
1178 clusters.

##  3 Methods

In this section, we describe the pipeline used to derive the cluster
demographics, motivated by the forward Bayesian modeling approach demonstrated
in Krumholz et al. (2019b). For reasons of brevity we only summarise the
method, and refer readers to Krumholz et al. for details.

###  3.1 Overview of the method

Given a set of unresolved photometric measurements of a cluster population,
how do we infer their underlying demographics? To answer this question, we
first propose a joint distribution of mass $M$, age $T$, and extinction
$A_{\mathrm{V}}$ for the cluster population, which we denote as
$f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})$, where $\boldsymbol{\theta}$
represents a vector of parameters describing the joint distribution $f$. For
example, if we assume the mass distribution of clusters is described by a
Schechter function $f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})\propto
M^{\alpha_{M}}\exp(-M/M_{\mathrm{break}})$, then $\boldsymbol{\theta}$
contains the slope $\alpha_{M}$ and break mass $M_{\mathrm{break}}$ of the
Schechter function. We describe the functional forms we consider for
$f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})$ and the parameters they involve
in Section 3.3.

We seek to derive the posterior distribution of $\boldsymbol{\theta}$. We do
so in the usual way for a Bayesian method, but writing the posterior as the
product of a prior and a likelihood function

$p_{\mathrm{post}}(\boldsymbol{\theta}\mid\\{\mathbf{m}\\})\propto\mathcal{L}(\\{\mathbf{m}\\}\mid\boldsymbol{\theta})p_{\mathrm{prior}}(\boldsymbol{\theta}).$
(1)

where $\mathcal{L}(\\{\mathbf{m}\\}\mid\boldsymbol{\theta})$ is the likelihood
for our set of $N_{\mathrm{obs}}$ photometric measurements. We defer a
discussion of priors to Section 3.3, and present details of our method of
calculating the likelihood function in Section 3.2, but to summarise the
latter here, we compute the likelihood function using a Gaussian mixture model
derived from a large library of synthetic star clusters; as we change
$\boldsymbol{\theta}$, we adjust the weights applied to this library, which
changes the value of $\mathcal{L}(\\{\mathbf{m}\\}\mid\boldsymbol{\theta})$.
Our method can therefore be summarised into the following steps.

  1. 1.

We generate a library composed of synthetic clusters with weights based on a
proposed distribution of cluster physical properties and the observational
completeness of the survey.

  2. 2.

Using the newly created library, we produce a synthetic distribution in
photometric space, which will be used to compare with the observations.

  3. 3.

We adjust the model parameters $\boldsymbol{\theta}$ to maximise the
resemblance between the observations and the photometric distribution, as
parameterised by the likelihood function. As we do so, we map out the
posterior distribution.

We carry out the final step of this procedure using a Markov Chain Monte Carlo
(MCMC) method as implemented in the software package emcee (Foreman-Mackey et
al., 2013). For all the calculations presented in this paper we use 4000
iterations of 100 walkers, discarding the first 300 iterations for burn-in;
both visual inspection of the posteriors and quantitative evaluation of the
auto-correlation time indicates that the chains are well-converged. We
describe the first two steps of the method – generating and re-weighting the
library, in Section 3.2.

###  3.2 Calculation of the likelihood function

As discussed above, we compute the likelihood function using the Gaussian
mixture model described in Krumholz et al. (2019b), which operates on a
library of $N_{\mathrm{lib}}$ synthetic clusters, each of which is
characterised by a mass $M$, age $T$, extinction $A_{\mathrm{V}}$, and a
vector of photometric magnitudes $\mathbf{m}$ in each of the filters used in
the observations. The full library is further described by a vector of
photometric bandwidths $\mathbf{h}$, which we set to 0.1 mag in all filters;
see Krumholz et al. (2015b) and Krumholz et al. (2019b) for detailed
discussion of the meaning of the bandwidth and the motivation for choosing
this value. We provide details regarding the library we use in this work in
Appendix A.

Given the library, Krumholz et al. (2019b) show that the likelihood function
for the parameters $\boldsymbol{\theta}$ can be written as (their equation 12)

$\mathcal{L}(\\{\mathbf{m}\\}\mid\boldsymbol{\theta})\propto\prod_{i=1}^{N_{\mathrm{obs}}}\left[\mathcal{A}(\boldsymbol{\theta})\sum_{j=1}^{N_{\mathrm{lib}}}w_{j}(\boldsymbol{\theta})\mathcal{N}(\mathbf{m}_{i}\mid\mathbf{m}_{j},\mathbf{h}^{\prime}_{i})\right],$
(2)

where $w_{j}(\boldsymbol{\theta})$ is the statistical weight of the $j$th
library cluster (which depends on the parameters $\boldsymbol{\theta}$ as
described below), $\mathbf{m}_{i}$ and $\mathbf{m}_{j}$ are the vectors of
magnitudes for the $i$th observed and $j$th library clusters, respectively,
$\mathcal{A}(\boldsymbol{\theta})=[\sum_{j=1}^{N_{\mathrm{lib}}}w_{j}(\boldsymbol{\theta})]^{-1}$
is a normalisation factor, and
$\mathcal{N}(\mathbf{x}\mid\mathbf{x}_{0},\boldsymbol{\sigma})$ is the
standard multidimensional Gaussian distribution centred at $\mathbf{x}_{0}$
and with width $\boldsymbol{\sigma}$, evaluated at position $\mathbf{x}$. The
quantity $\mathbf{h}_{i}^{\prime}$ is given by
$\mathbf{h}^{\prime}_{i}=(\mathbf{h}^{2}+\boldsymbol{\sigma}_{i}^{2})^{1/2}$,
where $\boldsymbol{\sigma}_{i}$ the observational error on the magnitude of
cluster $i$; the quantities $\mathbf{h}^{\prime}_{i}$, $\mathbf{h}$, and
$\boldsymbol{\sigma}_{i}$ are vectors with one element per filter, and the
expression for $\mathbf{h}^{\prime}_{i}$ should be understood as applying
separately to each filter. In practice we evaluate equation 2 numerically
using the cluster_slug module of the slug software suite, which implements a
tree-based order $N_{\mathrm{obs}}\ln N_{\mathrm{lib}}$ algorithm for
performing the calculation that is much faster than a naive brute force
evaluation, which would have a computational cost of order
$N_{\mathrm{obs}}N_{\mathrm{lib}}$

The quantity in square brackets in equation 2 is the distribution of
photometric magnitudes for the cluster library evaluated at the vector of
magnitudes $\mathbf{m}_{i}$ for the $i$th observed cluster. This distribution,
and thus the likelihood function as whole, depends on the parameters
describing the cluster population $\boldsymbol{\theta}$ only through the
weight functions $w_{j}(\boldsymbol{\theta})$ that describe the statistical
weight of each library cluster. The relationship between weights and
$\boldsymbol{\theta}$ is given by (equation 9 of Krumholz et al. 2019b)

$w_{j}(\boldsymbol{\theta})=P_{\mathrm{obs}}(\mathbf{m}_{j})\frac{f(M_{j},T_{j},A_{\mathrm{V},j}\mid\boldsymbol{\theta})}{p_{\mathrm{lib}}(M_{j},T_{j},A_{\mathrm{V},j})}.$
(3)

The denominator $p_{\mathrm{lib}(M_{j},T_{j},A_{\mathrm{V},j})}$ is the
distribution function describing the sampling density of the library (see
Appendix A), while the numerator
$f(M_{j},T_{j},A_{\mathrm{V},j}\mid\boldsymbol{\theta})$ is distribution of
the physical properties given $\boldsymbol{\theta}$, and the pre-factor
$P_{\mathrm{obs}}(\mathbf{m}_{j})$ is the probability that a cluster with
vector of magnitude $\mathrm{m}_{j}$ would be included in the LEGUS catalogue.
The terms in equation 3 can be intuitively understood as follows. The factor
$P_{\mathrm{obs}}(\mathbf{m}_{j})$ accounts for the fact that only a fraction
of clusters with magnitudes $\mathbf{m}_{j}$ will be observed due to the
completeness limits; we compute this probability as described in Section 2.3.
The denominator is the probability density for drawing clusters with a
particular combination of physical parameters $(M_{j},T_{j},A_{\mathrm{V},j})$
while constructing the library. Finally, the numerator represents the true
probability density for a cluster population described by the parameter set
$\boldsymbol{\theta}$, so that the ratio
$f(M_{j},T_{j},A_{\mathrm{V},j}\mid\boldsymbol{\theta})/p_{\mathrm{lib}}(M_{j},T_{j},A_{\mathrm{V},j})$
represents the factor by which we must up- or down-weight the library so that
clusters in the library have the same mass, age, and extinction distribution
as clusters in reality; if this weight factor is unity, then our library is
sampled from the same distribution of cluster properties as the real
population.

Intuitively, then, our method consists of iteratively adjusting the parameters
$\boldsymbol{\theta}$ and thus the weights $w_{j}(\boldsymbol{\theta})$ so as
to bring the predicted photometric distribution into as close agreement as
possible with the observed one. This will in turn adjust the mass and age
distributions, since these are determined by the same vector of parameters
$\boldsymbol{\theta}$ as the luminosity distribution. Again, we remind readers
that this is just an intuitive description of the underlying process; formal
proof that equation 2 is the correct likelihood function to accomplish this
adjustment, along with some details of how we handle complications like
clusters where some filters are missing due to the fields of view in the
different filters not being perfectly overlapping, is provided in Krumholz et
al. (2019b).

###  3.3 Demographic models and priors

Cluster demographics depend on cluster formation and destruction mechanisms,
and thus the set of candidate parametric model distributions
$f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})$ we consider is necessarily
informed by theoretical expectation. In this work, we consider the two most
prominent models: mass-independent destruction (MID; Section 3.3.1; Fall et
al. 2005, 2009; Chandar et al. 2015) and mass-dependent destruction/disruption
(MDD; Section 3.3.2; Lamers et al. 2005a; Gieles et al. 2007; Grudić et al.
2021). We couple both of these to a parametric model for the distribution of
extinctions (Section 3.3.3).

####  3.3.1 The mass-independent disruption (MID) model

In the MID model, the rate at which star clusters are destroyed is independent
of cluster mass. The cluster mass function therefore has the same shape at all
ages, and the distribution function
$f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})$ can be separated into two
distinct functions, one describing the mass distribution $p_{M}(M)$ and one
the age distribution $p_{T}(T)$. Given the observational evidence that the
mass function is a (possibly) truncated powerlaw, we will adopt a functional
form for $p_{M}(M)$ given by a Schechter function,

$p_{M}\left(M\right)\propto
M^{\alpha_{M}}\exp\left(-\frac{M}{M_{\mathrm{break}}}\right).$ (4)

For the age distribution, photometry cannot differentiate between a bound an
unbound cluster, so the age distribution is required to be flat for times that
are shorter than the physical time required for the stars in a cluster to
disperse; of course the distribution can also be flat out to older ages if the
mechanisms responsible for disruption do not begin until some time after a
cluster forms. Regardless of its physical, origin, we call the time at which
cluster disruption begins $T_{\mathrm{MID}}$. After this time, clusters will
disrupt, and we approximate the age distribution as a powerlaw. We therefore
have

$p_{T}(T)\propto\begin{cases}1,&\text{ if }T<T_{\text{MID }}\\\
\left(\frac{T}{T_{\text{MID }}}\right)^{\alpha_{T}},&\text{ if }T>T_{\text{MID
}}\end{cases}.$ (5)

Thus the joint mass-age distribution in the MID model is

$\frac{\mathrm{d}^{2}N}{\mathrm{~{}d}M\mathrm{~{}d}T}\propto
M^{\alpha_{M}}\exp\left(-\frac{M}{M_{\text{break
}}}\right)\max\left(T,T_{\text{MID }}\right)^{\alpha_{T}}.$ (6)

The MID model therefore has 4 free physical parameters that we place in our
vector $\boldsymbol{\theta}$: $\alpha_{M}$, $\log M_{\mathrm{break}}$,
$\alpha_{T}$, and $\log T_{\mathrm{MID}}$. We adopt flat priors on
$\alpha_{M}$ from $-4$ to 0 and $\alpha_{T}$ from $-3$ to 0, reflecting a
broad range around previous literature values; we will see that these choices
have little effect, as our MCMC never approaches these boundaries. The priors
on $\log M_{\mathrm{break}}$ and $\log T_{\mathrm{MID}}$ require somewhat more
thought. As for $\log M_{\mathrm{break}}$, we impose a flat prior from 2 to 7
because with we barely see clusters with masses more massive than
$10^{6.5}M_{\odot}$. The lower mass limit of 2 in log scale is to ensure a
reasonable MCMC walker range. For $\log(T_{\mathrm{MID}}/\mathrm{yr})$ we
impose flat priors from 5 to 10 based on physical plausibility. At the lower
end, clusters cannot disperse on less than a crossing timescale, and even the
densest clusters detected in LEGUS have crossing timescales well above
$10^{5}$ yr. The upper limit is roughly the age of the Universe.

####  3.3.2 The mass-dependent disruption (MDD) model

For the MDD model, clusters lose mass at a rate that varies as a pure power-
law function of their current mass (Gieles et al., 2007; Grudić et al., 2021),
$\mathrm{~{}d}M/\mathrm{~{}d}T\propto-M^{\gamma_{\mathrm{MDD}}}$. For such a
mass loss rate, the present-day mass of a cluster born with initial mass
$M_{i}$ at age $T$ is

$M=M_{\mathrm{i}}\left[1-\gamma_{\mathrm{MDD}}\left(\frac{M_{0}}{M_{\mathrm{i}}}\right)^{\gamma_{\text{MDD
}}}\frac{T}{T_{\mathrm{MDD},0}}\right]^{1/\gamma_{\text{MDD }}},$ (7)

where $T_{\mathrm{MDD},0}$ represents the time required for a cluster with
mass $M_{0}$ to have fully disrupted, defined as having reached a present-day
mass $M=0$. The joint mass-age distribution therefore obeys

$\frac{\mathrm{d}^{2}N}{\mathrm{~{}d}M\mathrm{~{}d}T}\propto\frac{\mathrm{d}^{2}N}{\mathrm{~{}d}M_{\mathrm{i}}\mathrm{~{}d}T}\frac{\mathrm{d}M}{\mathrm{~{}d}M_{\mathrm{i}}}.$
(8)

If the distribution of initial masses $M_{i}$ follows the Schechter function
form given by equation 4, then we have

$\frac{\mathrm{d}^{2}N}{\mathrm{~{}d}M\mathrm{~{}d}T}\propto
M^{\alpha_{M}}\eta^{\alpha_{M}+1-\gamma_{\text{MDD
}}}\exp\left(-\eta\frac{M}{M_{\text{break }}}\right),$ (9)

where

$\eta(M,T)=\left[1+\gamma_{\mathrm{MDD}}\left(\frac{M_{0}}{M}\right)^{\gamma_{\text{MDD
}}}\frac{T}{T_{\mathrm{MDD},0}}\right]^{1/\gamma_{\text{MDD }}}$ (10)

is the ratio of the initial and present-day cluster masses for a cluster of
present-day mass $M$ and age $T$.

The MDD model therefore has four free parameters: $\alpha_{M}$, $\log
M_{\mathrm{break}}$, $\gamma_{\mathrm{mdd}}$, and $T_{\mathrm{MDD,0}}$; note
that $M_{0}$ is not a separate parameter, because the mass-age distribution
depends only on the combination of parameters
$M_{0}^{\gamma_{\mathrm{MDD}}}/T_{\mathrm{MDD,0}}$. We therefore without loss
of generality choose $M_{0}=100$ M⊙ in what follows. However, for reader
convenience we will also report the commonly-used $t_{4}$ parameter, which is
simply the disruption time for a cluster of mass $10^{4}$ M⊙; this is given by
$t_{4}=10^{2\gamma_{\mathrm{MDD}}}T_{\mathrm{MDD,0}}$. We adopt the same
priors on $\alpha_{M}$ and $\log M_{\mathrm{break}}$ as in the MID model
(Section 3.3.1). For $\gamma_{\mathrm{MDD}}$, previous observational estimates
and $N$-body simulations give values in the range $0.6-0.7$ (Lamers et al.,
2010b), and we adopt broad priors that include this range: we take
$\gamma_{\mathrm{MDD}}$ to be flat from 0 to 1. These limits stem from
physical considerations: if $\gamma_{\mathrm{MDD}}<0$ then low-mass clusters
lose mass more slowly than massive ones, contrary to the physical expectations
of the model, while if $\gamma_{\mathrm{MDD}}\geq 1$ then no clusters ever
disrupt because there is no $T$ for which $M=0$. Finally, we adopt flat priors
on $\log(T_{\mathrm{MDD,0}}/\mathrm{yr})$ from 4 to 10; these limits are broad
enough not to matter, because none of our walkers ever approach them.

####  3.3.3 Dust extinction

The exact functional shape of the distribution of dust extinctions
$p_{A_{\mathrm{V}}}(A_{\mathrm{V}})$ is unknown, so we model it as non-
parametrically as possible. Following Krumholz et al. (2019b), we adopt a
simple piece-wise linear form over the range $A_{\mathrm{V}}=0-3$ mag
characterised by six free parameters ($p_{A_{\mathrm{V},i}},0\leq i\leq 6$) to
be fit, representing the value of the PDF at $A_{\mathrm{V}}$ = $i\,\Delta
A_{\mathrm{V}}$ mag with $\Delta A_{\mathrm{V}}=0.5$ mag. Thus the functional
form we adopt for the extinction distribution is

$\displaystyle p_{A_{\mathrm{V}}}(A_{\mathrm{V}})\propto$
$\displaystyle\begin{cases}p_{A_{\mathrm{V}},i}+\left(p_{A_{\mathrm{V}},i+1}-p_{A_{\mathrm{V}},i}\right)\left(\frac{A_{\mathrm{V}}}{\Delta
A_{\mathrm{V}}}-i\right),&i<\frac{A_{\mathrm{V}}}{\Delta A_{\mathrm{V}}}\leq
i+1\\\ 0,&A_{\mathrm{V}}>3\,\mathrm{mag}\end{cases}$

We treat $p_{A_{\mathrm{V},i}}$ for $i=0-5$ as parameters of our model to be
fit, with $p_{A_{\mathrm{V}},6}$ fixed by the requirement that the total area
under the PDF be unity. We set the priors on $p_{A_{\mathrm{V},i}}$ to be flat
for all values $>0$, subject to the requirement that $p_{A_{\mathrm{V},i}}>0$
remain positive for all $A_{\mathrm{V}}$.

Combining this with the mass and age distributions, our final functional form
to be fit is

$f(M,T,A_{\mathrm{V}}\mid\boldsymbol{\theta})\propto\frac{\mathrm{d}^{2}N}{\mathrm{~{}d}M\mathrm{~{}d}T}\,p_{A_{\mathrm{V}}}(A_{\mathrm{V}}),$
(12)

with $\mathrm{d}^{2}N/\mathrm{d}M\mathrm{~{}d}T$ given by equation 6 or
equation 9 for the MID or MDD models, respectively.

###  3.4 Model selection using Akaike weights

The method described thus far allows us to compute the posterior PDFs of the
model parameters $\boldsymbol{\theta}$ for both the MID and MDD models. To
determine whether the MID or MDD model provides a better fit and more accurate
description of the data, we use the Akaike information criterion (AIC) to
assess the fit quality provided by non-nested models. For models tested using
AIC, we first identify the walker with the highest likelihood, then compute,

$\mathrm{AIC}_{(\text{MID,MDD})}=2k-2\ln\hat{\mathcal{L}}_{\text{(MID,MDD)}},$
(13)

where $k$ represents the number of free parameters in the models, and
$\hat{\mathcal{L}}$ is the maximum of the likelihood function. In our case,
$k$ is 11 for both the MID and MDD models, with four parameters describing the
joint mass-age distribution ($\alpha_{M}$, $M_{\mathrm{break}}$, and either
$\alpha_{T}$ and $T_{\mathrm{MID}}$ or $\gamma_{\mathrm{MDD}}$ and
$T_{\mathrm{MDD,0}}$, six parameters representing the dust extinction shape,
and one extra parameter to describe the total number of clusters present in
the galaxy.

For the model comparison, we calculate the Akaike weights for either the MID
or MDD model as

$w_{\mathrm{MID/MDD}}=\frac{\frac{e^{-\Delta_{\mathrm{MID/MDD}}}}{2}}{\frac{e^{-\Delta_{\mathrm{MID}}}}{2}+\frac{e^{-\Delta_{\mathrm{MDD}}}}{2}}.$
(14)

The AIC measures the amount of information in the data preserved by a given
model., with the relative Akaike weight $w$ of one model indicating the
confidence level at which we can claim that it preserves more information than
the other models considered. The model with the highest Akaike weight is our
best fit and thus the preferred model.

##  4 Results

In this section, we report our fit parameters and model comparison results
(Section 4.1), along with comparisons between models and observed photometry
to verify that our best-fitting models do a reasonable job at reproducing the
observations (Section 4.2). To search for variations in cluster population
demographics with galactocentric radius, we also separately analyse clusters
in the inner and outer galaxy (Section 4.3).

Catalogue | Model | $w$ | $\alpha_{M}$ | $\log\left(M_{\mathrm{break}}/M_{\odot}\right)$ | $\alpha_{M_{4}}$ | $\alpha_{T}$ | $\log\left(T_{\mathrm{MID}}/\mathrm{yr}\right)$ | $\log(T_{\mathrm{MDD,0}}/\mathrm{yr})$ | $\gamma_{\mathrm{MDD}}$ | $t_{4}(\mathrm{Myr})$
---|---|---|---|---|---|---|---|---|---|---
All | MID | 0.9997 | $-1.03^{+0.10}_{-0.29}$ | $3.88^{+0.08}_{-0.06}$ | $-2.07^{+0.09}_{-0.28}$ | $-0.45^{+0.04}_{-0.07}$ | $6.78^{+0.42}_{-0.06}$ | $-$ | $-$ | $-$
MDD | 0.0003 | $-1.40^{+0.07}_{-0.08}$ | $4.00^{+0.05}_{-0.05}$ | $-$ | $-$ | $-$ | $6.65^{+0.10}_{-0.07}$ | $0.98^{+0.02}_{-0.04}$ | $268^{+37}_{-29}$
$\\{r\leq R_{50}\\}$ | MID | 1.0000 | $-0.54^{+0.18}_{-0.53}$ | $3.65^{+0.42}_{-0.07}$ | $-1.64^{+0.16}_{-0.45}$ | $-0.70^{+0.06}_{-0.19}$ | $6.78^{+0.29}_{-0.05}$ | $-$ | $-$ | $-$
MDD | < $10^{-11}$ | $-1.20^{+0.14}_{-0.13}$ | $3.87^{+0.08}_{-0.08}$ | $-$ | $-$ | $-$ | $6.42^{+0.11}_{-0.07}$ | $0.97^{+0.02}_{-0.05}$ | $151^{+22}_{-18}$
$\\{r>R_{50}\\}$ | MID | 0.9988 | $-1.12^{+0.13}_{-0.35}$ | $4.07^{+0.11}_{-0.08}$ | $-2.10^{+0.12}_{-0.32}$ | $-0.45^{+0.06}_{-1.41}$ | $6.85^{+1.44}_{-0.14}$ | $-$ | $-$ | $-$
MDD | 0.0012 | $-1.36^{+0.09}_{-0.09}$ | $4.13^{+0.07}_{-0.06}$ | $-$ | $-$ | $-$ | $6.75^{+0.20}_{-0.14}$ | $0.93^{+0.05}_{-0.10}$ | $274^{+52}_{-44}$
$\\{r\leq R_{c}\\}$ | MID | 1.0000 | $-0.90^{+0.09}_{-0.09}$ | $3.91^{+0.05}_{-0.05}$ | $-1.92^{+0.08}_{-0.07}$ | $-0.57^{+0.03}_{-0.04}$ | $6.76^{+0.05}_{-0.04}$ | $-$ | $-$ | $-$
MDD | < $10^{-11}$ | $-1.36^{+0.07}_{-0.07}$ | $4.05^{+0.04}_{-0.04}$ | $-$ | $-$ | $-$ | $6.48^{+0.07}_{-0.05}$ | $0.98^{+0.01}_{-0.03}$ | $183^{+19}_{-17}$
$\\{r>R_{c}\\}$ | MID | 0.9713 | $-1.47^{+0.33}_{-0.27}$ | $4.04^{+0.25}_{-0.17}$ | $-2.46^{+0.30}_{-0.22}$ | $-0.46^{+0.27}_{-1.74}$ | $8.23^{+0.42}_{-0.06}$ | $-$ | $-$ | $-$
MDD | 0.0287 | $-1.58^{+0.23}_{-0.23}$ | $4.06^{+0.31}_{-0.16}$ | $-$ | $-$ | $-$ | $8.58^{+0.88}_{-1.08}$ | $0.62^{+0.29}_{-0.43}$ | $2463^{+20082}_{-1905}$
$\\{p_{\mathrm{obs}}\geq 90\%\\}$ | MID | $-$ | $-1.42^{+0.11}_{-0.12}$ | $3.89^{+0.05}_{-0.5}$ | $-2.44^{+0.1}_{-0.1}$ | $-^{+0.03}_{-0.04}$ | $6.75^{+0.04}_{-0.04}$ | $-$ | $-$ | $-$
Table 1: Fitting results. The first column indicates the subset of the LEGUS
catalogue used in the fit, the second indicates the type of model fit (MID or
MDD), the third gives the Akaike weight for that model, and the remaining
columns give marginal posterior PDFs derived for each model parameter,
reported as ${q_{50}}^{+(q_{84}-q_{50})}_{-(q_{50}-q_{16})}$, where $q_{N}$
denotes the estimated $N^{\mathrm{th}}$ percentile. Thus our central values
are the median of the PDF, and the ranges shown correspond to the 68%
confidence interval. The co-rotation radius is denoted by $R_{c}$, and the
median galactocentric radius is denoted by $R_{50}$.

###  4.1 Full Catalogue Fits

We first analyse the full LEGUS NGC 628 catalogue using the method described
in Section 3. We summarise the marginal posteriors we derive on all model
parameters in the first block in Table 1. Our model comparison yields
$w(\mathrm{MID})=0.9997$ and $w(\mathrm{MDD})=3\times 10^{-4}$, suggesting
that the MID model does a substantially better job capturing the variations in
the data, and we will therefore focus on this as our preferred model from this
point forward.

Figure 2 shows the posterior PDFs of $\alpha_{M}$, $\log M_{\mathrm{break}}$,
$\alpha_{T}$, and $\log T_{\mathrm{MID}}$ we obtain from our fit. In this
plot, the extent of the axes reflects the full range of model parameters
allowed by our priors, so we can immediately read off where parameters are
well-constrained by the data, versus where they are unconstrained and occupy
the full range of values allowed by our priors. These figures have a few
noteworthy features. First, they suggest that cluster disruption in NGC 628 is
relatively mild; the most likely scenario is that disruption begins after
$T_{\mathrm{MID}}\approx 10^{7}$ yr but the is relatively slow
($\alpha_{T}\approx-0.5$), though an alternate scenario where cluster
disruption is faster but begins later ($\alpha_{T}\approx-1$ but
$T_{\mathrm{MID}}\approx 10^{8}$ yr) is marginally allowed. In the former
scenario, our median values of $\alpha_{T}\approx-0.45$ and
$\log(T_{\mathrm{MID}}/\mbox{yr})\approx 6.8$ correspond to cluster disruption
starting within $1-2$ crossing times of formation ($T_{\mathrm{MID}}\approx 6$
Myr), and thereafter occurring at a rate such that the cluster disruption time
is $\approx 2.2$ times the current cluster age (Krumholz et al., 2019a).
Similarly, the fits provide strong evidence that a truncated Schechter-form
CMF fits the data better than the pure power-law, with
$\log({M_{\text{break}}}/\mathrm{M}_{\odot})=3.88^{+0.08}_{-0.06}$ (68%
confidence interval), but that the powerlaw slope at lower masses is
relatively shallow, $\alpha_{M}\approx-1$. As with cluster disruption, there
is a tail of probability at higher $M_{\mathrm{break}}$ and steeper CMF
($\alpha_{M}\approx-2$), but this tail contains little probability mass.

Because our analysis method allows us to reach considerably smaller cluster
mass than many previous studies, and because we find a relatively small value
of $M_{\mathrm{break}}$, it is helpful to report an alternative measure of the
mass function slope that will be more readily comparable to literature
results. For this reason, we also compute the slope of the mass distribution
at a mass of $10^{4}$ M⊙,

$\alpha_{M_{4}}\equiv\left.\frac{d\log N}{d\log
M}\right|_{M=10^{4}M_{\odot}}=\alpha_{M}-\frac{10^{4}\,\mathrm{M}_{\odot}}{M_{\mathrm{break}}}.$

This is more directly comparable than $\alpha_{M}$ to the CMF slopes reported
in previous studies where the cluster sample is truncated at masses slightly
below $10^{4}$ M⊙ (e.g., Chandar et al., 2014; Adamo et al., 2017), and thus
the CMF slope is measured only at higher masses. Our results imply
$\alpha_{M_{4}}=-2.07_{-0.28}^{+0.09}$ (68% confidence), in good agreement
with the values of $\alpha_{M}\approx-2$ found in most previous studies (e.g.,
see Figure 5 of Krumholz et al. 2019a for a summary of observations).

Figure 2: Corner plot showing the one-dimensional and two-dimensional
histograms of the posterior PDFs of the parameters $\alpha_{M}$,
$\log{M_{\text{break}}}$, $\alpha_{T}$, and log $T_{\text{MID}}$ for the full
LEGUS NGC 628 catalogue. We omit the nuisance parameters describing the
distribution of extinctions $A_{V}$. The 1D blue histograms show one-
dimensional marginal PDFs for each parameter, while the colour maps show log
probability densities in various 2D cuts through parameter space. The
outermost contour level in the 2D plots is set so as to enclose 99% of the
samples, and black dots show individual MCMC samples outside this threshold.
The extent of each axis is rounded according to the prior range to achieve
optimal visualisation.

###  4.2 Photometric Comparisons

Before accepting the results of our fits, we must validate that our models
adequately recover the observed luminosity distribution, since matching this
distribution is the goal of our forward model. In order to visualise the
comparisons between the distributions of observations and the library, we
first plot the 1D cluster luminosity distributions in the five LEGUS bands,
showing both the measured distribution and the luminosity distribution we
predict using the model parameters. The model prediction follows immediately
from our expression for the likelihood function in terms of our Gaussian
mixture model (equation 2), and is simply

$p(m)\propto\sum_{j=1}^{N_{\mathrm{lib}}}w_{j}(\boldsymbol{\theta})\mathcal{N}(m\mid
m_{j},h),$ (15)

where $m$ is the magnitude in the filter of interest, $m_{j}$ is the magnitude
of the $j$th library cluster in that filter, and $h$ is the library bandwidth.
We show this comparison in Figure 3. As the figure shows, our best-fitting
model generally matches the observed luminosity distribution except for at the
most luminous end of the distribution.

Since our method attempts to fit the full five-dimensional photometric
distribution at once, we can also check for agreement between model and
observations in multiple dimensions. To this end, we present three colour-
colour plots illustrating the agreement between the observed and model-
predicted colour distributions in Figure 4, Figure 5 and Figure 6. As with the
luminosity distribution, we compute the colour distribution directly from the
Gaussian mixture model, using an expression analogous to equation 15. These
figures also show reasonably good agreement between the observed and model-
predicted colour distributions, and suggest that our best-fit model provides a
reasonable representation of the distribution of observations in 5D
photometric space.

The largest deviations between our model and the observations is for the most
luminous clusters at the bluest colours; we can see this from our slight
underprediction of the UV (F275W) and U luminosity functions in Figure 3, and
from the tail of observed clusters extending to higher UV - U (F275W-F336W)
than our models in Figure 4. One possible explanation for this discrepancy is
that our model library assumes Solar metallicity, while NGC 628 is slightly
sub-Solar (e.g., Berg et al., 2013, 2015). Krumholz et al. (2015b) show that
metallicity variations at this level have relatively little effect when
inferring cluster masses or ages, but it is conceivable that because our
library is somewhat more metal-rich than the data, it does not quite extend as
far to the blue as the observations. In any event, the effects of this
discrepancy on population-level statistics should be small, since they involve
only a handful of clusters out of our full catalogue of $>1000$.

Figure 3: Comparison between observed and model-predicted luminosity functions
in UV, U, B, V, I bands. Filled circles with error bars show the luminosity
distribution of the LEGUS catalogue, divided into 20 uniformly spaced bins;
distributions are normalised to have unit integral, and error bars show the
Poisson uncertainty in each bin. The black lines in the figure show model-
predicted photometric distributions computed using the $50^{\mathrm{th}}$
percentile values for all model parameters, while the grey bands around these
lines show the interquartile range for the models.  Figure 4: Comparison of
predicted and observed colour distributions for our best-fitting model. The
central panel shows the predicted colour distribution computed using the
$50^{\mathrm{th}}$ percentile values of all parameters as the background
heatmap, while blue circles show every $15^{\text{th}}$ cluster in the LEGUS
catalogue. The flanking panels compare 1D histograms of observed (red bars)
and predicted colour PDFs (blue bars). In order to visualize the Poisson error
of data, we plot the $84^{\text{th}}$ percentile upper limit on the
observational distribution as a solid red line.  Figure 5: Same as Figure 4,
except we show U-B versus B-V in colour-colour space here. Figure 6: Same as
Figure 4, except we show B-V versus V-I in colour-colour space.

###  4.3 Inner versus Outer Galaxy Clusters

Having presented the full catalogue fits, we now search for radial variations
in cluster demographics. For this purpose we assign every cluster in the LEGUS
catalogue a galactocentric radius, taking the centre of NGC 628 to be at
$\mathrm{RA}=24.17^{\circ}$, $\mathrm{DEC}=15.78^{\circ}$, and adopting an
inclination $i=25.2^{\circ}$ and position angle $\mathrm{PA}=25^{\circ}$ from
North (Grasha et al., 2015, 2017). We use the centre location, inclination,
and position angle to deproject the cluster coordinates from RA/DEC to a
galactocentric coordinate system. Once clusters have been assigned radii, we
divide the catalogue in two ways: first simply by making two sub-catalogues of
equal size containing the inner and outer halves of the sample, and second by
making two sub-catalogues containing clusters inside and outside the galactic
co-rotation radius. The advantage of the former approach is that it ensures
that we have equal statistical power in both regions, and thus do not miss
radial variations because our inner or outer galaxy sample is too small to see
them. The advantage of the latter approach is that co-rotation marks a
physically-motivated radius where we might expect to see a change in cluster
behaviour, as opposed to the galactocentric radius that marks an equal number
division, which is solely a function of the size and location of the LEGUS
pointings.

####  4.3.1 Equal Number Division

For our first analysis we divide the LEGUS cluster sample into two radial bins
each containing an equal number of clusters. The median radius that makes this
even division is $R_{50}=3.53$ kpc. We fit MID and MDD models to both the
subset of clusters with $r\leq R_{50}$ and those with $r>R_{50}$ and report
the results the second and third blocks of rows in Table 1. As with the full
catalogues, our model comparison strongly prefers the MID model in both the
inner and outer regions, and so we focus on it.

We compare the marginal posterior PDFs of $\alpha_{M_{4}}$ and
$M_{\mathrm{break}}$ in the inner and outer regions Figure 7. We notice that
the posterior mass distributions of the inner bin (indicated by the red
contours) are bimodal, corresponding to two probability peaks in the parameter
space, one at $\alpha_{M_{4}}\approx-1.8$ and
$\log(M_{\mathrm{break}}/\mathrm{M}_{\odot})\approx 3.8$, the second at
$\alpha_{M_{4}}\approx-2.5$ and
$\log(M_{\mathrm{break}}/\mathrm{M}_{\odot})\gtrsim 5$, while the posterior
for the outer region is closer to unimodal, centred at
$\alpha_{M_{4}}\approx-2$ and
$\log(M_{\mathrm{break}}/\mathrm{M}_{\odot})\approx 4.2$. Despite this
difference, however, the $2\sigma$ confidence contours for the two regions
overlap, so there is weak but not strongly convincing evidence for a
difference in mass function between the two regions.

Figure 8 provides an analogous plot for the two parameters – $\alpha_{T}$ and
$T_{\mathrm{MID}}$ – that describe the CAF. The qualitative conclusion to be
drawn from this Figure is similar to that for the mass PDFs as shown in Figure
7, i.e., the posterior PDFs for the inner and outer radial bins are largely
separated, but show a small region of overlap, in this case at
$T_{\mathrm{break}}\approx 6$ Myr and $\alpha_{T}\approx-0.7$. This result
suggests that we cannot rule out the possibility of two regions having
identical age distributions, but is suggestive of one. To the extent that
there is a difference, it is that the inner galaxy clusters show steeper
$\alpha_{T}$ and/or smaller $T_{\mathrm{MID}}$, suggesting the clusters closer
to the galactic centre undergo more severe disruption than those farther away.

####  4.3.2 Co-rotation Radius Division

An alternative way of separating the cluster population is into those inside
and outside the galactic co-rotation radius at $R_{c}=6.3$ kpc where the
spiral pattern and orbiting clusters move at equal speed (Sakhibov et al.,
2021). To the extent that encounters with molecular clouds or other structures
associated with the spiral arms influence cluster formation or disruption, we
might expect to see changes in cluster population demographics across co-
rotation. Making this division, we have 1159 clusters inside $R_{c}$ and 117
clusters outside $R_{c}$.

Figure 9 and Figure 10 show the inferred posteriors for the clusters located
inside and outside the co-rotation radius; we again concentrate on the MID
model, because model comparison prefers it to MDD. Here we see differences
between the inner and outer galaxy posterior PDFs that are qualitatively
consistent with those we observed for a division of the sample into two equal
parts. However, the small number of clusters in the outer galaxy sample when
we divide at co-rotation ensures that the posterior PDFs are very broad, and,
particularly for the age distribution, largely just reflect our priors. Given
the very broad outer galaxy posteriors, we cannot rule out the possibility
that the cluster demographics inside and outside co-rotation are the same.
Doing so would likely require significantly more outer-galaxy clusters than
the LEGUS catalogue provides.

Figure 7: 2D marginal posterior joint PDF of the parameters $\alpha_{M_{4}}$
and $\log M_{\mathrm{break}}$ that describe the CMF, for clusters inside (red)
and outside (blue) the median cluster galactocentric radius. Blue and red
contour lines in the central panel correspond to loci that enclose 68%, 95%,
and 99% of the MCMC samples, and thus correspond approximately to 1$\sigma$,
2$\sigma$ and 3$\sigma$ significance confidence levels. The histograms
flanking the central panel show the corresponding 1D marginal PDFs of
$\alpha_{M_{4}}$ and $\log M_{\mathrm{break}}$ for the inner (red) and outer
(blue) galaxy samples. Figure 8: Same as Figure 7, expect we show the marginal
PDFs of the parameters describing the CAF, $\alpha_{T}$ and $\log
T_{\mathrm{MID}}$, rather than the parameters describing the CMF. Figure 9:
Same as Figure 7, but for clusters divided into those inside and outside the
co-rotation radius $R_{c}$, rather than inside and outside the median radius
$R_{50}$.  Figure 10:  Same as Figure 8, but for clusters divided into those
inside and outside the co-rotation radius $R_{c}$, rather than inside and
outside the median radius $R_{50}$.  Figure 11: Same as Figure 7, but
comparing the results for the full LEGUS catalogue (blue) to those obtained
using only clusters in part of luminosity space where the estimated
completeness is $>90\%$. Figure 12: Same as Figure 11, expect we show age
distributions illustrated by 2D PDF of $\alpha_{T}$ versus
$\log\left(T_{\mathrm{MID}}/\mathrm{yr}\right)$.

##  5 Discussion

Here we discuss some implications of our findings. We begin with an additional
verification check on our method (Section 5.1), and then discuss our findings
regarding truncation of the cluster mass distribution (Section 5.2) and the
cluster age distribution (Section 5.3).

###  5.1 Method robustness against completeness uncertainty

In Section 2.3 we explain our method to derive the completeness function
$P_{\mathrm{obs}}(\mathbf{m})$ for the data and the library. Knowledge of this
function, coupled to our forward-modelling techniques, allows us to avoid the
traditional method of excluding clusters with masses $\lesssim 5000$ M⊙ and
ages $\gtrsim 200$ Myr, and in turn allowing us to use essentially all of the
LEGUS catalogue, rather than only $\approx 25\%$ of it as in previous work
(e.g., Adamo et al., 2017). However, our approach is sensitive to our
completeness estimates, which therefore represent a potential source of
systematic error. To test this possibility, we repeat the numerical
experiments with a more conservative approach of only including clusters where
we estimate the completeness is $\geq 90\%$. This reduces the sample to 1120
clusters, but ensures that $P_{\mathrm{obs}}(\mathbf{m})$ never varies by more
than 10%, and thus errors in it have little effect.

We compare the mass and age parameters we derive using this reduced catalogue
to those we obtain from the full catalogue in Figure 11 and Figure 12,
respectively. We find that the regions of allowed parameter space for both the
parameters describing mass ($\alpha_{M_{4}}$, $\log M_{\mathrm{break}}$) and
those describing age ($\alpha_{T}$, $\log T_{\mathrm{MID}}$) are consistent.
The probability peaks are slightly more concentrated for the high completeness
sample, but are almost entirely overlapping with the probability maxima we
obtain for the full catalogue; the results are consistent at the $1\sigma$
level. We see the similarity of the different physical parameters as shown in
Figures 11 and 12. Notably, we found consistent overlapping regions in these
three plots. In addition, the red contours fitted with high completeness
values have close probability peaks but are more concentrated than the entire
catalogue consisting of the fainter clusters. We therefore conclude that
possible systematic errors in the completeness function have at most a trivial
impact on the outcome.

###  5.2 CMF with a high mass truncation

Our results favour a CMF that is truncated at a relatively small mass,
$M_{\mathrm{break}}\approx 10^{4}$ M⊙, and are strongly incompatible with a
pure powerlaw with a slope $\alpha_{M}\approx-2$ – to the extent that our
results leave any room for a non-truncated mass function, they allow this
possibility only if the powerlaw slope $\alpha_{M}$ is significantly steeper
than $-2$. It is therefore interesting to compare our $M_{\mathrm{break}}$ to
values reported in the literature. Previous studies of CMF in spiral galaxies
such as M83 (Adamo et al., 2015), NGC 628 (Adamo et al., 2017), and M51 (Messa
et al., 2018) all found that Schechter functions provided better fits to the
CMF than pure powerlaws, but with exponential truncations in mass at $\approx
10^{5}M_{\odot}$, roughly an order of magnitude higher than our value. Indeed,
our break mass is roughly an order of magnitude smaller than the value
obtained by Adamo et al. (2017) from the same underlying photometric
catalogue.

One possible explanation for this difference is that the previous studies that
reported a higher $M_{\mathrm{break}}\sim 10^{5}M_{\odot}$ use mass-and age-
limited samples, $M>5000$ M⊙ and $T\leq 200$ Myr. These mass and age cuts
induce large uncertainties in the truncation mass estimate because they
require that the majority of the data be discarded. Indeed, it is worth
stressing that it would be essentially impossible for Adamo et al.’s study of
NGC 628 to have detected a break in the mass function at
$M_{\mathrm{break}}\approx 7600$ M⊙ as we find, simply because this value of
$M_{\mathrm{break}}$ is nearly the same as the mass at which they truncate
their sample; the could not have detected a change in CMF slope below this
mass. Consequently, CMFs reported in earlier studies should be understood as
describing only the CMF of the remaining heavy clusters after the bulk of the
population has been removed. This procedure may generate a high truncation
mass, whereas in our study these massive clusters, which are greatly
outnumbered by the low-mass ones, have much less weight in the fitting.
Consistent with this hypothesis, we note that the only two other published
studies of the extragalactic CMF that have mass sensitivity comparable to ours
(because they target galaxies close enough to resolve individual stars) also
report breaks in the mass function closer to $10^{4}$ than $10^{5}$ M⊙ –
Johnson et al. (2017) find $M_{\mathrm{break}}\approx 8500$ M⊙ in M31, while
Wainer et al. (2022) find $1.7\times 10^{4}$ M⊙ in M33. This strongly suggests
that in at least some cases there are real features in CMFs at masses below
the range accessible to traditional, non-Bayesian integrated light techniques.

Besides this statistical effect, the difference between our results and the
higher break masses found in some previous studies may also reflect real
changes in the cluster population due to galaxies’ physical conditions.
Johnson et al. (2017) propose a relation between the break mass and the star
formation rate density $\langle\Sigma_{\mathrm{SFR}}\rangle$, perhaps due to
both varying systematically with mean interstellar pressure; Wainer et al.
(2022) extend this proposal by adding new data. In order to place our fitted
break mass into this context, we add our value for $M_{\mathrm{break}}$ in NGC
628 to Wainer et al.’s compilation in Figure 13; for this purpose we adopt a
star formation rate density for NGC 628 from Santoro et al. (2022). As the
figure shows, our best-fit break mass is generally consistent with the
relation fit by Johnson et al.,

$\log\frac{M_{\mathrm{break}}}{\mathrm{M}_{\odot}}=(1.07\pm
0.10)\log\frac{\langle\Sigma_{\mathrm{SFR}}\rangle}{\mathrm{M}_{\odot}\;\mathrm{yr}^{-1}\;\mathrm{kpc}^{-2}}+(6.82\pm
0.20).$ (16)

This provides a possible physical explanation for why our break mass for NGC
628 is lower than the break masses typically found in galaxies like M51 and
M83: NGC 628 has a lower star formation rate per unit area, presumably
indicative of a lower interstellar pressure.

Figure 13: Comparison of break masses $M_{\mathrm{break}}$ and star formation
rate surface densities $\Sigma_{\mathrm{SFR}}$ for a sample of nearby
galaxies. Black points show the compilation of Wainer et al. (2022), while the
black line shows the best-fit linear relation obtained by Johnson et al.
(2017); they grey band shows the 84% confidence interval around the relation.
“M74(A2017)” denotes the truncated mass derived in Adamo et al. (2017) to make
it distinct from the mass cut-off reported in our study. The red point shows
our result for the full NGC 628 catalogue, with error bars indicating the 84%
confidence interval. We take our value of $\Sigma_{\mathrm{SFR}}$ for NGC 628
from Santoro et al. (2022).

###  5.3 Implications for cluster disruption

We find that mass-independent disruption (MID) is a better description of the
cluster age distribution than mass-dependent disruption (MDD) in all regions,
but that the disruption rate is relatively modest, $\alpha_{T}\approx-0.45$.
As with our results for the CMF, it is interesting to put this result into the
context of previous studies. As summarised by Krumholz et al. (2019a), the
results of these studies have tended to depend strongly on the method of
cluster catalogue construction, with authors who use “exclusive” catalogues
that exclude objects that do not have compact, round morphologies favouring
MDD and weak disruption, while those who use “inclusive” catalogues without
such morphological filters favouring MID and more rapid cluster disruption,
$\alpha_{T}\approx-1$. This difference arises because the morphological
filters tend to preferentially affect the youngest clusters, so applying them
flattens the age distribution. There is a great deal of debate in the
literature, which we shall not revisit here, about the extent to which
morphological filters artificially remove real star clusters from the sample,
versus the extent to which a failure to filter leads to catalogues including
large numbers of unbound associations that were never bound to begin with.

Because we include LEGUS class 3 sources in our catalogue, and thus are not
applying morphological filters, one would naively have expected our results to
more closely match those coming from the inclusive catalogues. However, our
actual result – MID but weak disruption – sits partway between the inclusive
and exclusive results, and this is likely at least in part because our method
allows us to investigate significantly older ages than in previous studies.
This makes us much less sensitive to the young ages where the uncertainty
about morphological filtering is most acute. Our statistical signal is not
dominated by the small number of clusters with ages $\lesssim 100$ Myr, but by
the much larger number of clusters at somewhat older ages, where
identification of clusters much less ambiguous since at these ages all bound
stellar systems are relaxed. Thus the situation here is somewhat analogous to
that for the CMF: it is not such much that our result is incompatible with
past work, as that our Bayesian forward-modelling approach provides access to
parts of parameter space that were not previously accessible.

It is also interesting to compare our results to those for NGC 628
specifically, for which Adamo et al. (2017) found mass-independent disruption
of clusters in the inner pointing NGC 628c, but mass-dependent disruption at
the outskirts in pointing NGC 628e. Since we choose to analyse the cluster
population as a whole or divide it by galactocentric radius, rather than by
pointing, the samples we analyse are not completely identical. However, since
the clusters belonging to NGC 628c (the inner region) constitute the great
majority of the population, our model selection result for NGC 628 as a whole
is consistent with the findings of Adamo et al.. There is possibly some
tension in the fact that we favour MID even when we divide our catalogues by
galactocentric radius, while Adamo et al. find MDD in NGC 628e, but since the
divisions are not exactly the same we cannot make this statement more
quantitatively. However, it does reinforce the point that our method has the
advantage that, since we are forward modelling the completeness in each field,
we can easily break up the clusters by any criterion of our choice to derive
the demographics in different parts of the galaxy.

##  6 Summary And Conclusion

We use a novel method described by Krumholz et al. (2019b) to derive the
demographics of the cluster population of the nearby spiral galaxy NGC 628
from unresolved photometry measured by the LEGUS survey (Calzetti et al.,
2015b). This method has the advantage of integrating diverse observational
data imaged using different filters or depths, without the need for binning or
imposing extreme completeness cuts. As a result, our method allows us to
analyze the roughly 1200 clusters in NGC 628 catalogued by the LEGUS survey,
as opposed to the past work that was restricted to include only the most
massive and youngest $\sim 30\%$ of the population. Our method also obviates
the need to assign a unique mass and age to every individual cluster of the
population, and therefore avoids the inevitable information loss incurred by
reducing complex, sometimes multiply-peaked posterior probability
distributions for individual cluster masses and ages to single values or even
Gaussians. We validate the method by comparing the photometric distributions
that emerge from our fits to the observed ones and finding good agreement,
indicating that our approach can identify models for the cluster population
that reproduce observations.

Our analysis results in three primary findings. First, we use the Akaike
information criterion to compare two classes of parametric models for the
time-evolution of the cluster population, mass-independent disruption (MID)
and mass-dependent disruption (MDD). We find that the MID model fits the data
better, both for the entire NGC 628 cluster catalogue and for all sub-
divisions of it that we consider. We thus conclude that cluster mass and age
distributions are at least approximately separable and that the disruption
time scales of the clusters have no dependence on the cluster masses.

Second, we find that cluster mass functions are well-described by Schechter
functions $dN/dM\propto M^{\alpha_{M}}\exp(-M/M_{\mathrm{break}})$, with a
well-constrained value of $M_{\mathrm{break}}\sim 10^{4}$ M⊙, indicating that
there is significant evidence for truncation in the cluster mass function.
Broadly, our finding indicates that the slope of the mass function changes
from $\alpha_{M}\sim-1$ at $10^{3}$ M⊙ to $\sim-2$ at $10^{4}$ M⊙ to much
steeper and nearly cut off at $\sim 10^{5}$ M⊙. Our truncation mass is
significantly lower than the one determined by Adamo et al. (2017) for NGC 628
but is in good agreement with those measured for M31 and M33 by Johnson et al.
(2017) and Wainer et al. (2022), respectively. These studies are based on
resolved stellar populations, and thus are similar to ours in that they can
probe significantly lower cluster masses than earlier studies such as Adamo et
al. (2017) based on integrated light that have found higher truncation masses.

Third, we find that the cluster age distribution in NGC 628 indicates that
cluster disruption begins when clusters are relatively young (age $\approx
5-15$ Myr), but that disruption is relatively mild thereafter. The age
distribution once disruption starts can be described by a simple power-law
with slope $\alpha_{T}\approx-0.4$. In past studies, the value of $\alpha_{T}$
has proven highly sensitive to the cluster catalogue type used, with studies
using an inclusive cluster catalogues (those that do not include filters based
on morphology) generally finding $\alpha_{T}\sim-1$ while those using
exclusive catalogues (which are filtered for morphology) finding
$\alpha_{T}\sim-0.3$. The significance of our study is that even though we are
using an inclusive catalogue, we still find a relatively shallow age slope.

The three conclusions above are derived from our analysis of all the clusters
in NGC 628 catalogued by LEGUS. However, we also group clusters by
galactocentric radius and search for differences in cluster demographics with
radius. We do not find strongly convincing evidence for any radial variations,
though there are suggestive hints that inner galaxy clusters are subject to
more severe disruption than those further out. Our ability to search for
radial variations is limited by the fairly small number of clusters available
at larger galactocentric radii, since the great majority of the available
sample comes from the LEGUS NGC 628c pointing, which targets the galactic
centre. Given the suggestive hints in our analysis, increasing the sample of
outer-galaxy clusters would be a worthwhile effort.

Looking forward, now that we have demonstrated the performance of our pipeline
on NGC 628, we are in a position to apply it to other cluster catalogues, both
those originating from LEGUS and from other studies. Analysis of a diverse
sample of galaxies will allow us to compare the demographics of cluster
populations in various galactic environments. This in turn opens up the
possibility of testing a range of theoretical models for cluster formation and
disruption and their variation with the galactic environment.

## Acknowledgements

JT acknowledges support from Kathryn Grasha’s Discovery Early Career
Researcher Award (DECRA) Fellowship (DE220100766). KG is supported by the
Australian Research Council through the Discovery Early Career Researcher
Award (DECRA) Fellowship DE220100766 funded by the Australian Government. MRK
acknowledges support from the ARC Future Fellowship and Laureate Fellowship
funding schemes, awards FT180100375 and FL220100020. This work is supported by
the Australian Research Council Centre of Excellence for All Sky Astrophysics
in 3 Dimensions (ASTRO 3D), through project number CE170100013. Based on
observations made with the NASA/ESA Hubble Space Telescope, obtained at the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
These observations are associated with program #13364.

## Data Availability

The cluster catalog data underlying this article are publicly available online
at MAST. The software suite, tools, and LEGUS catalogues required to perform
the MCMC optimization, as well as the output chains, are publicly accessible
on Bitbucket.

## References

  * Adamo et al. (2015) Adamo A., Kruijssen J. M. D., Bastian N., Silva-Villa E., Ryon J., 2015, MNRAS, 452, 246
  * Adamo et al. (2017) Adamo A., et al., 2017, ApJ, 841, 131
  * Bastian et al. (2005) Bastian N., Hempel M., Kissler-Patig M., Homeier N. L., Trancho G., 2005, A&A, 435, 65
  * Bastian et al. (2012) Bastian N., et al., 2012, MNRAS, 419, 2606
  * Baumgardt & Kroupa (2007) Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589
  * Berg et al. (2013) Berg D. A., Skillman E. D., Garnett D. R., Croxall K. V., Marble A. R., Smith J. D., Gordon K., Kennicutt Jr. R. C., 2013, ApJ, 775, 128
  * Berg et al. (2015) Berg D. A., Skillman E. D., Croxall K. V., Pogge R. W., Moustakas J., Johnson-Groh M., 2015, ApJ, 806, 16
  * Bialopetravičius et al. (2019) Bialopetravičius J., Narbutis D., Vansevičius V., 2019, A&A, 621, A103
  * Calzetti et al. (2015a) Calzetti D., et al., 2015a, AJ, 149, 51
  * Calzetti et al. (2015b) Calzetti D., et al., 2015b, AJ, 149, 51
  * Chandar et al. (2014) Chandar R., Whitmore B. C., Calzetti D., O’Connell R., 2014, ApJ, 787, 17
  * Chandar et al. (2015) Chandar R., Fall S. M., Whitmore B. C., 2015, ApJ, 810, 1
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Dalcanton et al. (2012) Dalcanton J. J., et al., 2012, ApJS, 200, 18
  * Dinnbier & Walch (2020) Dinnbier F., Walch S., 2020, MNRAS, 499, 748
  * El-Badry et al. (2017) El-Badry K., Weisz D. R., Quataert E., 2017, MNRAS, 468, 319
  * Fall & Chandar (2012) Fall S. M., Chandar R., 2012, ApJ, 752, 96
  * Fall & Zhang (2001) Fall S. M., Zhang Q., 2001, ApJ, 561, 751
  * Fall et al. (2005) Fall S. M., Chandar R., Whitmore B. C., 2005, ApJ, 631, L133
  * Fall et al. (2009) Fall S. M., Chandar R., Whitmore B. C., 2009, ApJ, 704, 453
  * Fensch et al. (2019) Fensch J., et al., 2019, A&A, 628, A60
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fouesneau et al. (2014) Fouesneau M., et al., 2014, ApJ, 786, 117
  * Gieles & Renaud (2016) Gieles M., Renaud F., 2016, MNRAS, 463, L103
  * Gieles et al. (2006) Gieles M., Portegies Zwart S. F., Baumgardt H., Athanassoula E., Lamers H. J. G. L. M., Sipior M., Leenaarts J., 2006, MNRAS, 371, 793
  * Gieles et al. (2007) Gieles M., Lamers H. J. G. L. M., Portegies Zwart S. F., 2007, ApJ, 668, 268
  * Grasha et al. (2015) Grasha K., et al., 2015, ApJ, 815, 93
  * Grasha et al. (2017) Grasha K., et al., 2017, ApJ, 842, 25
  * Grudić et al. (2021) Grudić M. Y., Kruijssen J. M. D., Faucher-Giguère C.-A., Hopkins P. F., Ma X., Quataert E., Boylan-Kolchin M., 2021, MNRAS, 506, 3239
  * Gusev et al. (2014) Gusev A. S., Egorov O. V., Sakhibov F., 2014, MNRAS, 437, 1337
  * Hills (1980) Hills J. G., 1980, ApJ, 235, 986
  * Howard et al. (2017) Howard C. S., Pudritz R. E., Harris W. E., 2017, MNRAS, 470, 3346
  * Johnson et al. (2017) Johnson L. C., et al., 2017, ApJ, 839, 78
  * Kruijssen et al. (2012) Kruijssen J. M. D., Pelupessy F. I., Lamers H. J. G. L. M., Portegies Zwart S. F., Bastian N., Icke V., 2012, MNRAS, 421, 1927
  * Krumholz et al. (2015a) Krumholz M. R., Fumagalli M., da Silva R. L., Rendahl T., Parra J., 2015a, MNRAS, 452, 1447
  * Krumholz et al. (2015b) Krumholz M. R., et al., 2015b, ApJ, 812, 147
  * Krumholz et al. (2019a) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019a, ARA&A, 57, 227
  * Krumholz et al. (2019b) Krumholz M. R., Adamo A., Fumagalli M., Calzetti D., 2019b, MNRAS, 482, 3550
  * Lamers et al. (2005a) Lamers H. J. G. L. M., Gieles M., Portegies Zwart S. F., 2005a, A&A, 429, 173
  * Lamers et al. (2005b) Lamers H. J. G. L. M., Gieles M., Bastian N., Baumgardt H., Kharchenko N. V., Portegies Zwart S., 2005b, A&A, 441, 117
  * Lamers et al. (2010a) Lamers H. J. G. L. M., Baumgardt H., Gieles M., 2010a, MNRAS, 409, 305
  * Lamers et al. (2010b) Lamers H. J. G. L. M., Baumgardt H., Gieles M., 2010b, MNRAS, 409, 305
  * Larsen (2002) Larsen S. S., 2002, AJ, 124, 1393
  * Lee et al. (2022) Lee J. C., et al., 2022, ApJS, 258, 10
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Linden et al. (2017) Linden S. T., et al., 2017, ApJ, 843, 91
  * Maschberger & Kroupa (2009) Maschberger T., Kroupa P., 2009, MNRAS, 395, 931
  * Messa et al. (2018) Messa M., et al., 2018, MNRAS, 473, 996
  * Mok et al. (2020) Mok A., Chandar R., Fall S. M., 2020, ApJ, 893, 135
  * Olivares E. et al. (2010) Olivares E. F., et al., 2010, ApJ, 715, 833
  * Piskunov et al. (2009) Piskunov A. E., Kharchenko N. V., Schilbach E., Röser S., Scholz R. D., Zinnecker H., 2009, A&A, 507, L5
  * Sakhibov et al. (2021) Sakhibov F., Gusev A. S., Hemmerich C., 2021, MNRAS, 508, 912
  * Santoro et al. (2022) Santoro F., et al., 2022, A&A, 658, A188
  * Wainer et al. (2022) Wainer T. M., et al., 2022, ApJ, 928, 15
  * Webb et al. (2019) Webb J. J., Reina-Campos M., Kruijssen J. M. D., 2019, MNRAS, 486, 5879
  * da Silva et al. (2012) da Silva R. L., Fumagalli M., Krumholz M., 2012, ApJ, 745, 145
  * de Meulenaer et al. (2014) de Meulenaer P., Narbutis D., Mineikis T., Vansevičius V., 2014, A&A, 569, A4

##  Appendix A The cluster library

The Krumholz et al. (2019b) approach described in the main text requires a
library of synthetic clusters. We construct a library of $4\times 10^{7}$
clusters for this purpose using the slug stochastic stellar population
synthesis code (da Silva et al., 2012; Krumholz et al., 2015a). Construction
of the library involves two steps. In the first, we choose a mass $M$, age
$T$, and extinction $A_{\mathrm{V}}$ for each cluster from a specified library
distribution $p_{\mathrm{lib}}(M,T,A_{\mathrm{V}})$. We choose this library
distribution to be

$p_{\mathrm{lib}}\left(M,T,A_{\mathrm{V}}\right)\propto
p(M)p(T)p\left(A_{\mathrm{V}}\right),$ (17)

with

$\displaystyle p(M)$ $\displaystyle\propto$
$\displaystyle\begin{cases}(M/10^{2}\,\mathrm{M}_{\odot})^{-1},\quad&10^{2}<M/\mathrm{M}_{\odot}\leq
10^{5}\\\
10^{-3}(M/10^{5}\,\mathrm{M}_{\odot})^{-2},&10^{5}<M/\mathrm{M}_{\odot}\leq
10^{7}\\\ 0,&\text{otherwise}\end{cases}$ $\displaystyle p(T)$
$\displaystyle\propto$
$\displaystyle\begin{cases}T^{-1},\quad&10^{5}<T/\mbox{yr}<1.5\times
10^{10}\\\ 0,&\text{otherwise}\end{cases}$ $\displaystyle
p\left(A_{\mathrm{V}}\right)$ $\displaystyle\propto$
$\displaystyle\begin{cases}\text{const},\quad&0<A_{\mathrm{V}}<3\text{ mag}\\\
0,&\text{otherwise}\end{cases}$

This distribution maximises sampling of clusters with low masses and young
ages, where stochastic effects are strongest and thus the largest number of
samples is needed to fully characterise the luminosity distribution.

The second step is to generate photometric magnitudes for each cluster. To do
so, we run slug using the MIST v1.0 (Choi et al., 2016) tracks at Solar
metallicity, coupled to the “starburst99” treatment of stellar atmospheres
(Leitherer et al., 1999). We draw stars from a Chabrier IMF using slug’s “stop
nearest” sampling policy. Full explanations of these choices are provided in
Krumholz et al. (2015a), and in the slug documentation available at
https://slug2.readthedocs.io/. However, Krumholz et al. (2015b) show that,
within reasonable limits, making alternative choices of IMF, metallicity, and
stellar tracks leads to relatively modest systematic variations in the cluster
properties that slug infers.


# Emergent Dynamics in Neural Cellular Automata

Yitao Xu, Ehsan Pajouheshgar, Sabine Süsstrunk

School of Computer and Communication Sciences, EPFL, Switzerland
{yitao.xu, ehsan.pajouheshgar, sabine.susstrunk<EMAIL_ADDRESS>
###### Abstract

Neural Cellular Automata (NCA) models are trainable variations of traditional
Cellular Automata (CA). Emergent motion in the patterns created by NCA has
been successfully applied to synthesize dynamic textures. However, the
conditions required for an NCA to display dynamic patterns remain unexplored.
Here, we investigate the relationship between the NCA architecture and the
emergent dynamics of the trained models. Specifically, we vary the number of
channels in the cell state and the number of hidden neurons in the MultiLayer
Perceptron (MLP), and draw a relationship between the combination of these two
variables and the motion strength between successive frames. Our analysis
reveals that the disparity and proportionality between these two variables
have a strong correlation with the emergent dynamics in the NCA output. We
thus propose a design principle for creating dynamic NCA.

## Introduction

Cellular Automata’s (CA) ability to create complex patterns has fascinated
researchers for many years (Von Neumann et al.,, 1966; Gardner,, 1970; Wolfram
and Gad-el Hak,, 2003). Neural Cellular Automata (NCA), a more recent
development (Mordvintsev et al.,, 2020), introduces a trainable variant of
conventional CA. The emergence of complex and life-like behaviors, such as
spontaneous motion and self-regeneration from simple local update rules, makes
NCA a compelling candidate for studying artificial life. NCAs have been
applied to a variety of problems, including modeling morphogenesis
(Mordvintsev et al.,, 2020; Sudhakaran et al.,, 2021, 2022), texture synthesis
(Niklasson et al., 2021b, ; Pajouheshgar et al., 2023b, ; Pajouheshgar et al.,
2023a, ), and collective modeling of intelligent behavior (Randazzo et al.,,
2020; Nadizar et al.,, 2022), demonstrating the effectiveness of NCAs in
modeling complex biological and natural processes.

Recent studies (Pajouheshgar et al., 2023b, ; Pajouheshgar et al., 2023a, )
have particularly focused on utilizing emergent NCA dynamics to create moving
patterns. Pajouheshgar et al., 2023b  highlight that in NCA models trained for
texture synthesis, the collective of cells exhibits spontaneously moving
patterns. They propose a training schema to give structure to this random
emergent motion and create dynamic textures. However, understanding the
necessary conditions that give rise to such a spontaneous motion remains
unexplored. Here, we investigate the relationship between emergent NCA
dynamics and its hyperparameters, including the number of channels in the cell
states ${C}$, and the number of hidden neurons $D$, in the MLP. We provide
qualitative and quantitative results to show that NCA loses its dynamics as
the disparity or ratio of $D$ and $C$ decreases. Based on our findings, we
propose an NCA design principle to retain the dynamic properties.

## Preliminaries

### Neural Cellular Automata

In NCA, cells exist on a 2D grid, where each cell’s state is represented by a
$C$ dimensional vector. We refer to $C$ as the Channel. Given a 2D grid of
size $H\times W$, all cells collectively define the state of the system
$\mathbf{S}\in\mathbb{R}^{C\times H\times W}$. An NCA iteratively updates cell
states $\mathbf{S}$ over time by alternating between Perception and Adaptation
stages. A single update step is visualized in Figure 1.

Figure 1: A single update step of NCA.

In the Perception stage, cells perceive their local 9-point Moore neighborhood
via fixed convolution kernels. The typical choice (Mordvintsev et al.,, 2020;
Niklasson et al., 2021b, ; Pajouheshgar et al., 2023b, ) of the convolution
kernels are Sobel filters in the $x$ and $y$ directions, $K_{x},K_{y}$, and a
Laplacian filter $K_{lap}$. Moreover, the cell itself is also counted as part
of its neighborhood through convolution using an identity filter $K_{id}$. The
output of the Perception stage is then fed into an MLP, denoted $F_{\theta}$,
where $\theta$ denotes the set of trainable parameters in the Adaptation
stage. $D$, also known as Hidden, is the number of hidden neurons in the MLP.
Finally, the output of the MLP is multiplied by a stochastic binary mask,
$M\sim Bernoulli(0.5)$, to ensure asynchronous cell update, as this property
is important for NCAs to behave well (Niklasson et al., 2021a, ). In summary,
the NCA update rule can be expressed as:

$\mathbf{S}^{t+1}=\mathbf{S}^{t}+F_{\theta}(\mathbf{S}^{t},(K_{x},K_{y},K_{lap})*\mathbf{S}^{t})\odot
M$ (1)

### Motion Quantification

The patterns created by the consecutive application of the NCA update rule
display emergent and spontaneous motion. To quantify the strength of this
motion, we utilize the pre-trained Optical Flow estimation network from
(Tesfaldet et al.,, 2018). Given two consecutive frames from the sequence of
patterns synthesized by an NCA, let $O_{ij}\in\mathbb{R}^{2}$ be the estimated
optical flow at the pixel $i,j$. We define the motion strength $\Psi$ as the
norm of optical flow averaged across all pixels:

$\Psi=\frac{1}{HW}\sum^{H}_{i=1}\sum^{W}_{j=1}||O_{ij}||_{2}.$ (2)

Figure 2: $Y$ axis shows the average motion strength $(\Psi)$. (a): $X$ axis
shows Hidden - Channel $(D-C)$. Correlation significance: $p<0.002$. Small $Y$
values are excluded in the correlation fitting for a better visualization.
With all values: $r=0.6225,p<0.0001$. (b): $X$ axis shows Hidden / Channel
$(\frac{D}{C})$. Correlation significance: $p<0.0001$. Figure 3: Visualization
of frames during the test of NCA using two different numbers of channels $C$.
Top row: $C=24,D=96$. Bottom row: $C=96,D=128$. $O(S^{t},S^{t+1})$: the optic
flow (Baker et al.,, 2011) between the two frames. $|S^{t}-S^{t+1}|$: the
absolute difference between the two frames. Better viewed in video format.

## Experiments and Results

For different values of $C$ in `range(8,128,8)` , and $D$ in
`range(16,128,16)` , we train NCAs on four different textures in the DTD
dataset (Cimpoi et al.,, 2014): bubbly_0101, chequered_0121, interlaced_0172,
cracked_0085, resulting in a total of 512 trained models. All models are
trained for 6000 epochs to ensure convergence. We average the motion strength
$\Psi$ across 100 frames for each NCA configuration. Each frame corresponds to
$T$ NCA steps. We find that $T=32$ creates sufficiently close yet distinct
enough frames for accurate optic flow prediction. Moreover, we observe that
all NCAs become stable after around 100 frames. Therefore, we first let the
NCA generate 100 frames and then perform the measurement on the second 100
frames. Figure 2-(a) shows the relationship between motion strength and $D-C$
for all NCA configurations. We observe that once $D<C$, the emergent motion
completely disappears and as the difference $D-C$ increases, the NCA becomes
more dynamic. However, there are still some NCA models with $D>C$ that are not
dynamic. We further investigate this problem by plotting the ratio instead of
the difference between $D$ and $C$ in Figure 2-(b). We notice that NCA barely
has motion when $\frac{D}{C}<2$, which accounts for the case when there is no
motion, but $D>C$. We visualize the output of two NCA models, A: $C=24,D=96$
and B: $C=96,D=128$, in Figure 3 to demonstrate an example of dynamic and non-
dynamic NCA configurations. Although NCA B satisfies $D>C$, the ratio between
$D$ and $C$ is too small for the model to permit dynamic patterns.

## Discussion and Conclusion

We demonstrate that the number of hidden neurons should be larger than the
number of channels and that their ratio should be at least greater than 2.0
for an NCA to display emergent dynamics. Our experiments effectively disprove
the statement that “NCA can underfit or overfit to the static target image and
thus become non-dynamic”, since both small and large NCAs can show no motion
when $\frac{D}{C}<2$. Without the dynamic properties, NCA is prone to get
stuck in local minima where artifacts appear in the generated texture, as
shown in Figure 3. Moreover, the lack of dynamics will render NCA no longer
suitable for certain applications, such as dynamic texture synthesis and
biological pattern simulation. Therefore, in practice, if we favor dynamic
NCAs without artifacts, it is critical to ensure the Hidden $D$ is larger by a
factor of two or more than the Channel $C$.

## References

  * Baker et al., (2011) Baker, S., Scharstein, D., Lewis, J., Roth, S., Black, M. J., and Szeliski, R. (2011).  A database and evaluation methodology for optical flow.  International Journal of Computer Vision, 92(1):1–31.
  * Cimpoi et al., (2014) Cimpoi, M., Maji, S., Kokkinos, I., Mohamed, S., , and Vedaldi, A. (2014).  Describing textures in the wild.  In Proceedings of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR).
  * Gardner, (1970) Gardner, M. (1970).  Mathematical games.  Scientific american, 222(6):132–140.
  * Mordvintsev et al., (2020) Mordvintsev, A., Randazzo, E., Niklasson, E., and Levin, M. (2020).  Growing neural cellular automata.  Distill.  https://distill.pub/2020/growing-ca.
  * Nadizar et al., (2022) Nadizar, G., Medvet, E., Nichele, S., and Pontes-Filho, S. (2022).  Collective control of modular soft robots via embodied spiking neural cellular automata.  arXiv preprint arXiv:2204.02099.
  * (6) Niklasson, E., Mordvintsev, A., and Randazzo, E. (2021a).  Asynchronicity in neural cellular automata.  In Artificial Life Conference Proceedings 33, volume 2021, page 116. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info ….
  * (7) Niklasson, E., Mordvintsev, A., Randazzo, E., and Levin, M. (2021b).  Self-organising textures.  Distill, 6(2):e00027–003.
  * (8) Pajouheshgar, E., Xu, Y., Mordvintsev, A., Niklasson, E., Zhang, T., and Süsstrunk, S. (2023a).  Mesh neural cellular automata.  arXiv preprint arXiv:2311.02820.
  * (9) Pajouheshgar, E., Xu, Y., Zhang, T., and Süsstrunk, S. (2023b).  Dynca: Real-time dynamic texture synthesis using neural cellular automata.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 20742–20751.
  * Randazzo et al., (2020) Randazzo, E., Mordvintsev, A., Niklasson, E., Levin, M., and Greydanus, S. (2020).  Self-classifying mnist digits.  Distill.  https://distill.pub/2020/selforg/mnist.
  * Sudhakaran et al., (2021) Sudhakaran, S., Grbic, D., Li, S., Katona, A., Najarro, E., Glanois, C., and Risi, S. (2021).  Growing 3d artefacts and functional machines with neural cellular automata.  In Artificial Life Conference Proceedings 33, volume 2021, page 108. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info ….
  * Sudhakaran et al., (2022) Sudhakaran, S., Najarro, E., and Risi, S. (2022).  Goal-guided neural cellular automata: Learning to control self-organising systems.  arXiv preprint arXiv:2205.06806.
  * Tesfaldet et al., (2018) Tesfaldet, M., Brubaker, M. A., and Derpanis, K. G. (2018).  Two-stream convolutional networks for dynamic texture synthesis.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 6703–6712.
  * Von Neumann et al., (1966) Von Neumann, J., Burks, A. W., et al. (1966).  Theory of self-reproducing automata.  IEEE Transactions on Neural Networks, 5(1):3–14.
  * Wolfram and Gad-el Hak, (2003) Wolfram, S. and Gad-el Hak, M. (2003).  A new kind of science.  Appl. Mech. Rev., 56(2):B18–B19.