Measuring Star-Formation Histories, Distances, and Metallicities with Pixel Color-Magnitude Diagrams I: Model Definition and Mock Tests

TL;DR

This paper introduces an advanced GPU-accelerated pixel color-magnitude diagram technique for accurately measuring star formation histories, metallicities, distances, and dust properties of galaxies, even in semi-resolved regimes, with high precision.

Contribution

The paper presents a new Python package, pcmdpy, with improved physical modeling and computational speed, enabling detailed stellar population analysis of galaxies up to 100 Mpc away.

Findings

The method can recover SFH, metallicity, distance, and dust extinction from mock data.

GPU acceleration makes the analysis roughly 7 times faster.

The technique is effective for galaxies within 10 Mpc and applicable up to 100 Mpc.

Abstract

We present a comprehensive study of the applications of the pixel color-magnitude diagram (pCMD) technique for measuring star formation histories (SFHs) and other stellar population parameters of galaxies, and demonstrate that the technique can also constrain distances. SFHs have previously been measured through either the modeling of resolved-star CMDs or of integrated-light SEDs, yet neither approach can easily be applied to galaxies in the "semi-resolved regime". The pCMD technique has previously been shown to have the potential to measure stellar populations and star formation histories in semi-resolved galaxies. Here we present Pixel Color-Magnitude Diagrams with Python (pcmdpy), a GPU-accelerated package that makes significant computational improvements to the original code and including more realistic physical models. These advances include the simultaneous fitting of distance, modeling a Gaussian metallicity-distribution function, and an observationally-motivated dust model. GPU-acceleration allows these more realistic models to be fit roughly 7x faster than the simpler models in the original code. We present results from a suite of mock tests, showing that with proper model assumptions, the code can simultaneously recover SFH, [Fe/H], distance, and dust extinction. Our results suggest the code, applied to observations with HST-like resolution, should constrain these properties with high precision within 10 Mpc and can be applied to systems out to as far as 100 Mpc. pCMDs open a new window to studying the stellar populations of many galaxies that cannot be readily studied through other means.

Tables

Table 1. Table 1: Modeling a pCMD requires an assumed physical model for each of: metallicity, star-formation history, dust extinction, and distance. The physical model implementations in PCMDPy and their parameters are listed here. Free-parameters are fit to data, while hyper-parameters are set prior to fitting.

No.	Model Name	Free Parameters	Hyper-parameters
Metallicity Models
M1	Single $\left[\mathrm{Fe}/H\right]$	$\left[\mathrm{Fe}/H\right]$
M2	Gaussian MDF (Fixed-Width)	$\left[\mathrm{Fe}/H\right]$	${\sigma }_{\left[\mathrm{Fe}/H\right]}$
M3	Gaussian MDF	$\left[\mathrm{Fe}/H\right]$, ${\sigma }_{\left[\mathrm{Fe}/H\right]}$
Star Formation History Models
S1	Single Stellar Population	$\log$$N_{\mathrm{pix}}$, $\log \mathrm{age}$
S2	Constant SFR	$\log$$N_{\mathrm{pix}}$
S3	Tau Model	$\log$$N_{\mathrm{pix}}$, $\tau$
S4	Delayed-Tau Model	$\log$$N_{\mathrm{pix}}$, $\tau$
S5	Non-Parametric SFH	$\log {\mathrm{SFH}}_0,\log {\mathrm{SFH}}_1,\cdots$	$N_{\mathrm{SFH}}$, Bin Edges
Dust Extinction Models
E1	Fixed Dust Screen	$\log E\left(B-V\right)$	$F_{\mathrm{dust}}$
E2	Log-Normal Dust Screen (Fixed-Width)	$\log E\left(B-V\right)$	$F_{\mathrm{dust}}$, ${\sigma }_{\mathrm{dust}}$
E3	Log-Normal Dust Screen	$\log E\left(B-V\right)$, ${\sigma }_{\mathrm{dust}}$	$F_{\mathrm{dust}}$
Distance Models
D1	Fixed Distance		${\mu }_d$
D2	Variable Distance	${\mu }_d$

Table 2. Table 2: Global hyper-parameters of the PCMDPy modeling procedure. These parameters remain fixed throughout fitting. Does not include parameters of the sampling algorithm (see Dynesty documentation).

Parameter Name (Default)	Num Params	Description
Filters	$N_{\mathrm{bands}}$	Specify which observational filters to simulate. Includes zero-points and PSF models.
$N_{\mathrm{im}}$ (512)	1	Size of simulated image-plane ($N_{\mathrm{im}}$ pix $\times$ $N_{\mathrm{im}}$ pix)
$N_{\mathrm{age}}$(21)	1	Number of isochrones to draw from in a SFH
IMF (Salpeter)	1	The stellar initial-mass function
Sky Level (0)	$N_{\mathrm{bands}}$	Level of background sky noise to add to each band
Exposure Time	$N_{\mathrm{bands}}$	Exposure time of each image (important for modeling shot noise)
Downsampling (5)	1	Factor to downsample isochrone points
Hess Binning (0.05, 0.05)	2	Width of Hess bins used to compute Log-Likelihood

Equations

$${\mathrm{N}}_{\mathrm{pix}}$$

$$\mathcal{L}$$

$${\sigma }_i$$

$$\mathcal{C}$$

Full text

galaxies: stellar content; galaxies: photometry; techniques: photometric;

††software: This research has made use of NASA’s Astrophysics Data System, as well as the following software packages: PyCUDA (Klöckner et al., 2012), Dynesty (Speagle in prep), NumPy (van der Walt et al., 2011), Matplotlib (Hunter, 2007), IPython (Pérez & Granger, 2007), Jupyter (Kluyver et al., 2016), SciPy (Jones et al., 2001), Pandas (McKinney, 2010), and Astropy (The Astropy Collaboration et al., 2013, 2018).

1 Introduction

The mass build-up of galaxies is sensitive to many physical processes. Mergers, AGN activity, supernovae feedback, and the accretion of pristine gas from cosmological filaments will all affect a galaxy’s star-formation history (SFH) and chemical enrichment, and it is a major goal of modern astrophysics to constrain their relative impacts. Hydrodynamical simulations (e.g., Hopkins et al., 2014; Vogelsberger et al., 2014) are one powerful way to study the evolution of galaxies, but their results require reliable observational constraints on the diversity of SFHs and stellar populations observed in the universe.

To date, measurements of stellar populations and SFHs largely rely on two distinct techniques, each with limitations on the systems they can be applied to and the robustness of the results. These approaches can be considered in a single framework through considering the typical number of stars per resolution element, $\mathrm{N_{pix}}$ (van Dokkum & Conroy, 2014; Conroy & van Dokkum, 2016).

The first technique, resolved star photometry, typically requires $\mathrm{N_{pix}}\lesssim\mathcal{O}(10^{-1})$ in order to fully resolve each individual star in a system. The fluxes of these stars in at least 2 bands are converted into a color-magnitude diagram (CMD), and the stellar populations are derived from fits to stellar evolution models (e.g., Dolphin, 2002; Weisz et al., 2011; Lewis et al., 2015; Williams et al., 2015). As long as stars down to the oldest main-sequence turnoff are resolved, this method recovers highly precise and robust measurements of the SFH, and is considered the gold standard.

In practice, the number of systems where the oldest main-sequence turnoff can be resolved is small, even with the resolution of HST’s Advanced Camera for Surveys (ACS). Several dozen dwarf galaxies in the local group have sufficiently low stellar densities to allow well-measured SFHs (e.g., Weisz et al., 2011, 2014). But even our nearest massive neighbor, M31, is so crowded that only rare, massive main-sequence stars (Lewis et al., 2015) can be individually resolved in most fields, limiting SFH recovery to relatively recent star formation. In the inner regions of M31, with $\mathrm{N_{pix}}\sim\mathcal{O}(10^{1})$, SFH can be measured using the resolved red giant branch or red clump stars, but the results are subject to significant systematics due to disagreement between models of these phases of stellar evolution (Williams et al., 2017), with the oldest ages the most uncertain. The inner bulge of M31 is so crowded that not even RGB stars can be resolved with HST, making any SFH measurement impossible with this technique.

The second well-established method, SED modeling, analyzes the integrated light of all stars in a (typically entirely unresolved) galaxy. The broadband photometry and/or spectra are modeled using stellar population synthesis techniques to recover stellar populations and SFHs (e.g., Walcher et al., 2011; Conroy, 2013). These methods typically assume a fully populated initial mass function (IMF), ignoring Poisson fluctuations of rare, evolved stars. Such a regime requires $\mathrm{N_{pix}}\gtrsim\mathcal{O}(10^{6})$, usually a valid approximation for distant, unresolved galaxies.

Yet even if such fluctuations can be ignored, the light from the oldest, lowest-mass main-sequence stars is usually dwarfed by the light of rare, evolved stars, an effect known as "outshining", which can lead to an underestimate of the total stellar mass and oldest ages of star formation (Maraston et al., 2010; Pforr et al., 2012; Sorba & Sawicki, 2015). The SFHs recovered from SED modeling are also highly dependent on the underlying stellar evolution models (e.g., Conroy, 2013; Hunt et al., 2018), with systematic uncertainties much larger than for resolved star SFHs. Carnall et al. (2018) and Leja et al. (2018) also study the prior bias introduced from assuming various functional forms for the SFH.

In between the realm of resolved stars and integrated light, there is a substantial volume of the universe in the so-called semi-resolved regime (van Dokkum & Conroy, 2014; Conroy & van Dokkum, 2016). This regime can be defined loosely as $\mathcal{O}(10^{1})\lesssim\mathrm{N_{pix}}\lesssim\mathcal{O}(10^{6})$, where stars cannot be individually resolved due to crowding, but surface-brightness fluctuations due to Poisson sampling of rare, bright stars in each pixel cannot be ignored. This semi-resolved regime extends from $\sim$1 Mpc out to nearly 100 Mpc for massive galaxies observed with HST.

The surface brightness fluctuations (SBF; Tonry & Schneider, 1988) distance technique is one example of a method for studying galaxies in the semi-resolved regime. In this approach, distances to a galaxy are estimated by measuring the scale of the pixel-to-pixel surface-brightness fluctuations, studied in Fourier space in order to isolate the instrumental PSF. The underlying stellar populations are largely treated as nuisance parameters and accounted for through calibration. Other examples include the "disintegrated" light analysis of stellar halos (Mould, 2012), fluctuation spectroscopy (van Dokkum & Conroy, 2014), and pixel-level time variability (Conroy et al., 2015).

A third technique for measuring SFHs and stellar populations, specifically applicable to semi-resolved galaxies, was introduced in Conroy & van Dokkum (2016). The pixel color-magnitude diagram (pCMD) method measures the fluxes in multiple (typically two) photometric bands within every pixel in an image, plotting the resulting magnitudes in a color-magnitude diagram. The surface-brightness fluctuations across the image represent real variations in the number of rare, bright stars at each pixel. Therefore, the distribution of these pixel fluxes (the pCMD) holds important information about the underlying stellar populations. Critically, as we show in this work, pCMD distributions show distinct sensitivities to many key physical parameters, including metallicity, SFH, distance, and dust content. With a proper accounting of important observational artifacts, including PSF convolution, sky subtraction, and shot noise, these key parameters can be inferred through a forward-modeling procedure.

The pCMD technique can therefore be considered a generalization of SBF, as in this work we demonstrate that both distances and stellar populations can be simultaneously measured using the fluctuations in multiple bands of well-aligned photometry. SBF is still an important rung on the distance ladder to this day, used to derive distances to nearby dwarf galaxies (Cohen et al., 2018; Carlsten et al., 2019) and massive galaxies alike (Greene et al., 2019). SBF also remains the best distance estimator to the host galaxy of GW170817, the first binary neutron star merger detected in gravitational waves (Cantiello et al., 2018). Modeling the pCMDs of systems such as these offers the prospect of constraining (rather than assuming) their stellar populations, while more robustly accounting for their uncertainties in deriving distances.

This work presents the first results from Pixel Color-Magnitude Diagrams with Python (PCMDPy), an open-source package developed to fully implement the pCMD fitting method and allow for modeling pCMDs with more realistic, flexible physical models. The package extends significantly on the original framework outlined in Conroy & van Dokkum (2016), allowing for a variety of metallicity, dust extinction, and star-formation history models, and for the first time extends the pCMD method to simultaneously fit for galactic distances. The package is written in Python and hosted on GitHub111https://github.com/bacook17/pcmdpy, final API and documentation still in the development stages, and posterior estimation is implemented with the dynamic nested sampling code dynesty (Speagle in prep)222https://github.com/joshspeagle/dynesty. The computationally challenging simulation procedure has been accelerated using graphics processing units (GPUs), allowing for more accurate and complex models to be generated more rapidly and for fits to converge in less time.

This paper is outlined as follows. In Section 2 we provide an outline of the PCMDPy code. In Section 3 we describe a suite of mock tests where simulated pCMDs generated with the code are then fit to study the code’s constraining power, and conclude the paper in Section 4.

2 Pixel Color-Magnitude Diagrams with Python (PCMDPy)

We describe here the procedure for generating and modeling a pCMD. Section 2.1 discusses the general approach and primary assumptions made, while Section 2.2 details the physical models implemented in PCMDPy, including metallicity, star-formation history, distance, and dust extinction. Section 2.3 demonstrates each free parameter in these models has a unique effect on pCMD distributions. In Section 2.4 we describe the computational architecture of the code. Section 2.5 describes our likelihood model for comparing pCMDs and the nested sampling approach for fitting model posteriors over these free parameters.

2.1 Overview of the pCMD Technique

A pixel color-magnitude diagram (pCMD, Conroy & van Dokkum, 2016) represents the distribution in pixel-by-pixel photometry of a galaxy or galactic region, projected into magnitude space. A pCMD is generated by measuring the flux in each pixel of two bands of photometric images and converting those fluxes to apparent magnitudes, representing their distribution in color-magnitude space as a Hess diagram. When doing so, we are necessarily discarding any spatial correlations between the pixels. This contrasts with the SBF method, where the Fourier transform preserves spatial information and can therefore isolate the effect of the PSF. The distribution of the pCMD is thus shaped both by the physical properties (stellar populations) of the stars residing in the image and by the dust extinction and the optical properties of the telescope, such as the PSF and photometric noise.

The pCMD technique endeavours to constrain the underlying stellar populations through a forward modeling approach: given a model for stellar photometry and our knowledge of the telescope’s optical properties, we can create a simulated pCMD derived from a specified set of populations and compare it to real data. This procedure is summarized visually in Figure 1. Comparing simulated pCMDs over a range of stellar populations allows us to infer the most likely populations to have generated in the data, as long as our stellar and observational models are reasonable approximations of the truth.

In practice, the possible combinations of stars residing in a given image is limitless, requiring us to make simplifying assumptions for the problem to be tractable. We adopt a probabilistic approach: we assume there is a global, underlying distribution of stars within the image, from which each pixel represents a unique random realization. Therefore, we care only about the overall distribution of stellar populations and pixel magnitudes (the pCMD), not about matching each specific pixel to a model of the stars it contains.

Specifically, we assume that the number of stars in a pixel is drawn from a Poisson distribution, with mean number of stars constant across the image and equal to the free-parameter $\mathrm{N_{pix}}$. Therefore, the surface-brightness fluctuations across pixels are a true signal, representing the Poisson noise in stars per pixel, the magnitude of which is determined by $\mathrm{N_{pix}}$ and the stellar populations. $\mathrm{N_{pix}}$ can be estimated from the column density $\mathrm{N}_{\mathrm{col}}$ of stars, the observed spatial resolution $\theta$, and the distance $\mathrm{d}$ to the source, and is approximately:

[TABLE]

The representative values above assume a galaxy with average stellar density of 1 star per cubic parsec and a thickness of 1 kpc ($\mathrm{N}_{\mathrm{col}}=10^{3}\mathrm{pc}^{-2}$), and observations with the spatial resolution of HST-ACS.

To determine the properties of those $\approx\mathrm{N_{pix}}$ stars in each pixel, we assume a model for the distribution of their ages (the star-formation history model, SFH), metallicities (given as iron abundance, $[\mathrm{Fe}/\mathrm{H}]$), and initial masses (an initial mass function, IMF). This is shown in Panel A of Figure 1.

In our approach, we divide the age-$[\mathrm{Fe}/\mathrm{H}]$-mass space into discrete bins, assigning each point a weight according to the models above and normalizing the weights to equal $\mathrm{N_{pix}}$. For every pixel in the simulated image (of size $\mathrm{N}_{\mathrm{im}}\times\mathrm{N}_{\mathrm{im}}$ pixels), we randomly draw the number of stars in each bin according to a Poisson distribution with the given weights333Given many Poisson-distributed numbers $n_{i}$ with weight parameters $\lambda_{i}$, the sum $N=\sum_{i}n_{i}$ is also Poisson-distributed, with weight parameter $\lambda=\sum_{i}\lambda_{i}$.. We then derive the absolute magnitudes of each star in the observed filters from stellar evolution models (isochrones), although this could also be implemented with well-calibrated observational catalogues. This is shown in Panel B of Figure 1.

The intrinsic fluxes of each star are attenuated by dust extinction, potentially by both foreground (Milky Way) and source (host galaxy) dust, requiring reddening curves and an assumed model for the dust abundance. The stellar magnitudes are converted to fluxes according to the distance to the host galaxy and are summed into their respective pixels. These steps are Panels C and D of Figure 1.

Finally, the telescope and instrumental signatures must be accounted for (Panel E of Figure 1). The raw images are convolved with models for the point-spread function (PSF) in each filter444The PSF convolution breaks the assumption that each pixel represents an independent draw from an underlying distribution, as the fluxes in neighboring pixels are highly correlated. This makes writing down an exact likelihood or computing it using probabilistic programming intractable.. Sky flux can be added at this stage. Poisson shot noise is added, assuming the data are in electrons, such that noise $=\sqrt{\mathrm{counts}}$. The resulting pCMD is computed by converting the pixel counts to magnitudes (Panel F of Figure 1).

We additionally note that it is possible to both measure and simulate pCMDs in any arbitrary number of photometric bands. In this work we only consider the case of two filters, primarily due to the challenge of comparing pCMD distributions in more than 2 dimensions.

2.2 PCMDPy Model Choices

The PCMDPy physical model for a region contains four primary model components: a metallicity model, a SFH model, a dust extinction model, and a distance model. When fitting a pCMD, the free parameters of the fit all correspond to one of these four components, and we assume flat priors over each. Each component (with the exception of distance) could be modeled in a variety of complex ways, and we detail here the implementations available in PCMDPy. The functional forms included here may not necessarily be adequate representations of the true shapes of the distributions of interest in particular galaxies, but they are designed to provide reasonable approximations. In Section 3.4, we discuss the effects of incorrect model assumptions.

All of these physical model components are summarized in Table 1. All other parameters or constants assumed in the model are considered hyper-parameters, and are held fixed throughout the fit. These are summarized in Table 2.

Three metallicity models are available, which specify the metallicity-distribution function (MDF, in terms of $[\mathrm{Fe}/\mathrm{H}]$) of the stars in each pixel.

1. M1.

Single [Fe/H]: All stars have the same metallicity, equal to the single free parameter $[\mathrm{Fe}/\mathrm{H}]$. 2. M2.

Gaussian MDF (Fixed-Width): Stellar metallicities are drawn from a Gaussian distribution, with one free parameter corresponding to the mean ($[\mathrm{Fe}/\mathrm{H}]$). The standard-deviation ($\sigma_{[\mathrm{Fe}/\mathrm{H}]}$) is a hyper-parameter and is held fixed. 3. M3.

Gaussian MDF: As in M2, but $\sigma_{[\mathrm{Fe}/\mathrm{H}]}$ is a second free parameter.

We adopt the non-rotating isochrones of the MIST project, v1.2 (Choi et al., 2016). The metallicities provided in the MIST isochrones are discretized to a grid with spacing 0.25 dex. For model M1, we interpolate the isochrones between these grid points to recover metallicities not on the grid. For models M2 and M3, we use only the metallicities on the grid, and weight their abundances by the specified Gaussian distribution. Grid points with weights less than $1\%$ are removed, resulting in 5 to 10 metallicity points for most Gaussian models.

Five star-formation history models are implemented, which specify the distribution of stars as a function of age. The number of stars per pixel $\mathrm{N_{pix}}$ is either an explicit free-parameter of the model, or computed as the sum of the SFH.

1. S1.

Single Stellar Population: All stars are of the same age. The two free parameters are the age of stars (in log years) and $\mathrm{N_{pix}}$. 2. S2.

Constant SFR: Assumes constant star-formation rate (SFR) at all times. The one free parameter is $\mathrm{N_{pix}}$, to which the total SFH sums. 3. S3.

Tau SFH: An exponentially-falling SFR ($\mathrm{SFR}\propto e^{-t/\tau}$). The two free parameters are $\tau$ (in Gyr) and $\mathrm{N_{pix}}$. 4. S4.

Delayed-Tau SFH: A linearly-rising SFR followed by an exponential falloff ($\mathrm{SFR}\propto te^{-t/\tau}$). The two free parameters are $\tau$ (in Gyr) and $\mathrm{N_{pix}}$. 5. S5.

Non-Parametric SFH: The star-formation history is binned into several (by default: 5) independent bins, within which the SFR is constant. The free parameters are $\{\log\mathrm{SFH}_{i}\}$, the logarithm of total star-formation in each bin, in units of stars-per-pixel. The number and edges of the bins are hyper-parameters.

With the exception of model S1, stellar ages are discretized using a grid of ages. By default, the grid uses 21 equally-spaced bins in log age from 1 Myr to 14 Gyr, but this is a hyper-parameter $\mathrm{N_{age}}$ that can be adjusted. The age points are taken as the midpoints of each bin ($10^{6.1}$, $10^{6.3}$, $\cdots$, $10^{10.1}$ years), and the star-formation rate is assumed to be constant within the bin.

When simulating complex models represented by many metallicity and age points (ex: M2+S3), we downsample the MIST isochrones in mass by a factor of $5$ (a hyper-parameter). This downsampling factor improves computation time dramatically while not significantly affecting the resulting pCMDs, as confirmed through internal tests.

In the current implementation of PCMDPy, SFH and metal abundance are modeled independently: the metallicity of a star does not depend on its age. In reality, the abundances of stars are known to evolve with age as previous generations enrich the ISM out of which stars form. Future work could model both jointly, but this is outside the scope of the current work.

We adopt a Salpeter IMF (Salpeter, 1955) by default, but a Kroupa IMF (Kroupa & Pavel, 2001) is also implemented in PCMDPy. In addition to determining the distribution of stars by mass, the IMF also determines the conversion from $\mathrm{N_{pix}}$ to $\mathrm{M_{pix}}$, the stellar mass formed in each pixel. The total stellar mass in the image, $\mathrm{M_{\star}}$, requires a solution for SFH to account for stellar mass loss from post-MS stellar evolution.

Three dust extinction models are implemented, which determine the dust extinction in units of the reddening parameter, $\log\mathrm{E}(\mathrm{B}-\mathrm{V}){}$. We adopt the framework of Dalcanton et al. (2015), who studied of the dust extinction in M31. See their Section 3 for details.

1. E1.

Constant Dust Screen: All pixels have a constant amount of extinction, equal to the one free parameter: $\log\mathrm{E}(\mathrm{B}-\mathrm{V})$. 2. E2.

Log-Normal Dust Screen (Fixed-Width): The dust extinction in each pixel is drawn from a log-normal distribution, with one free parameter: the median extinction $\log\mathrm{E}(\mathrm{B}-\mathrm{V})$. The dimensionless width-parameter $\sigma_{\mathrm{dust}}$ is a hyper-parameter and is held fixed. 3. E3.

Log-Normal Dust Screen: As in E2, except $\sigma_{\mathrm{dust}}$ is an additional free parameter.

Each dust model assumes a single, thin screen of dust. The geometry of the screen is specified by a hyper-parameter $\mathrm{F}_{\mathrm{dust}}$, which determines the fraction of stars that are reddened by dust. $\mathrm{F}_{\mathrm{dust}}=1.0$ represents a foreground screen of dust (all stars are reddened), while our default choice of $\mathrm{F}_{\mathrm{dust}}=0.5$ corresponds to a mid-plane disk of dust, with half the stars reddened and half unobscured. This model assumes that in practical applications, any foreground dust from the Milky Way can be accounted for in data reduction. Future work could extend to more complex dust geometries, as preliminary tests indicate neither of these models may be sufficient to model the complex and dense dust lane structures found in disk galaxies at the scale of interest. We convert $\log\mathrm{E}(\mathrm{B}-\mathrm{V}){}$ to magnitudes of extinction in each optical band using the $R_{V}=3.1$ reddening law from Schlafly & Finkbeiner (2011, their table 6).

The final physical model component is distance. We implement two distance models, simply representing whether or not distance (in units of distance modulus, $\mu_{d}$) is included in the fit:

1. D1.

Distance Fixed: Distance is assumed known. $\mu_{d}$ is a fixed hyper-parameter. 2. D2.

Distance Free: $\mu_{d}$ is a free parameter.

A physical model is specified with one of each of these model components. We will occasionally refer to the entire model used to simulate a pCMD using the alpha-numerical codes above. For example, our "fiducial-$\tau$" model (Section 3.1) is described as M2+S3+E2+D2. This translates to a Gaussian MDF (with fixed $\sigma_{[\mathrm{Fe}/\mathrm{H}]}$), a Tau SFH, a Log-Normal dust screen (with fixed $\sigma_{\mathrm{dust}}$), and distance free.

A major hyper-parameter of the fitting technique is $\mathrm{N}_{\mathrm{im}}$, the 1D size of the simulated image plane (the total number of pixels in the simulated pCMD is therefore $\mathrm{N}_{\mathrm{im}}^{2}$). $\mathrm{N}_{\mathrm{im}}$ does not have to match the size of the data being fit, as we can compare the relative number of pixels in the pCMD. Larger simulated images are more computationally expensive but reduce the inherent stochasticity of the likelihoods (see Section 2.5) by providing more samples of the rare fluctuations in surface-brightness. We choose a default image size of $\mathrm{N}_{\mathrm{im}}=512$, which we find to be the optimal size for reducing stochasticity while allowing for convergent fits in a reasonable time.

We must match the observational conditions under which the data were taken as closely as possible, requiring several well-calibrated hyper-parameters representing each observed filter. This includes the exposure time in each filter, the photometric zero-points, and a model for the PSF. We have specifically modeled observations with the HST-ACS camera, but the model can be generally applied to ground-based observations as well. The exposure time is taken from the FITS header of the imaging data. Zero-points are computed using the PySynphot (Lim et al., 2015) package, as detailed in the ACS Handbook555http://www.stsci.edu/hst/acs/analysis/zeropoints.

PSF models are taken from the Tiny Tim (Krist et al., 2011) web interface666http://www.stsci.edu/hst/observatory/focus/TinyTim. To simulate sub-pixel PSF effects, we subdivide each image into a $4\times{}4$ grid and apply to each a different PSF convolution, shifted by fractions of a pixel. Tests showed our approach is a reasonable approximation to the effects of truly sub-sampling each pixel and applying a sub-pixel PSF model, and is substantially faster.

2.3 Sensitivity of pCMDs to Model Parameters

In Figure 2, we show that the detailed structure of pCMDs are sensitive to each of the physical parameters of interest, in ways that allow constraints on each from only two bands of semi-resolved photometry. In each panel, we show a simulated pCMD for a "baseline model" in grey, with a comparison for a model with one physical parameter changed superimposed in black. The baseline model has a single-$[\mathrm{Fe}/\mathrm{H}]$ metallicity model (M1, $[\mathrm{Fe}/\mathrm{H}]=-0.5$), a Tau SFH (S3, $\mathrm{N_{pix}}=10^{3},\;\tau=3$ Gyr), constant dust screen (D1, $\log\mathrm{E}(\mathrm{B}-\mathrm{V}){}=-0.5$), and a distance of 1 Mpc ($\mu_{d}=25$).

Of particular interest are the effects of $\mathrm{N_{pix}}$ (proxy for stellar surface density) and distance. The upper-left column of Figure 2 shows that in addition to changing the average luminosity, $\mathrm{N_{pix}}$ also affects the dispersion of the pCMD distribution, due to the increase in Poisson surface-brightness fluctuations. We compare this to variations in distance modulus in the upper-right column, which simply shifts the distribution in the vertical (luminosity) direction, with no effect on color. We show variations in $\mathrm{N_{pix}}$ and $\mu_{d}$ that each result in the same average luminosity; the dispersion of the distribution can be used to break the degeneracy between luminosity and distance. This represents a reformulation of the SBF method: as the same physical system is moved towards larger distance, surface-brightness remains constant but the magnitude of fluctuations decreases as $\mathrm{N_{pix}}$ rises. This hints at the utility of pCMDs to simultaneously recover distances and stellar populations for galaxies in the semi-resolved regime. We demonstrate this conclusively in Section 3.

Varying metallicity shifts the peak of the distribution and has a notable effect on the slope of the upper-right wing. The pixels in that region contain red-giant branch (RGB) stars, and the pCMD is therefore sensitive to the metallicity-dependent slope of the RGB (Choi et al., 2016). The effects of increasing dust broadens the overall distribution but leaves the slope of the RGB-feature relatively intact. As we show in Section 3.2, the model is able to constrain dust and $[\mathrm{Fe}/\mathrm{H}]$, although there remains a slight degeneracy between the two.

Changing the age parameter, $\tau$, has significant effects on the blue-wing of the pCMD. This represents the relative abundance of young, massive main-sequence stars. These hot stars are only abundant when there has been recent star-formation (high $\tau$), while only the old main-sequence and RGB stars contribute when recent star-formation is suppressed (low $\tau$).

2.4 Computational Infrastructure

Even when downsampling the isochrones by a factor of 5, simulating a pCMD with the simplest SSP model requires $\mathcal{O}(10^{8})$ random Poisson calls for a $512\times 512$ image. Increasing the complexity and realism of a model by incorporating more age and metallicity points increases the computation time nearly linearly with the number of isochrones. With 21 age points and $\sim 5$ metallicity points, the more realistic physical models described in Section 2.2 require $\mathcal{O}(10^{10})$ Poisson draws, and quickly become computationally infeasible on a traditional CPU.

PCMDPy includes an optional GPU-accelerated backend, which dramatically reduces the computational time required to simulate a pCMD. Given the isochrones representing a particular physical model, each GPU thread independently samples the stars from those isochrones into an individual pixel. This accelerates the simulation time of an individual model pCMD by a factor of $\sim 30\times$ compared to the CPU implementation. Figure 3 shows the computation time required to simulate a pCMD of a single stellar population (model M1+S1) as a function of $\mathrm{N}_{\mathrm{im}}$. The CPU tests were run on an Intel Xeon E5-2620 processor (2.1 GHz) and the GPU tests on an Nvidia Tesla K20Xm777These are the default GPUs available to us. The latest generation of Nvidia GPUs (the Tesla V100) could lead to an additional $\sim 4\times$ speed-up..

The original code of Conroy & van Dokkum (2016) simulated an $\mathrm{N}_{\mathrm{im}}=256$ pCMD for a relatively simple physical model (M1+S2+E1+D1) in $\approx 1$s. The same computation in PCMDPy takes only $\approx 0.25$s on an Nvidia Tesla K20Xm. For our preferred model with a Gaussian MDF and 21 age bins (see Section 3.1), PCMDPy simulates a pCMD with $\mathrm{N}_{\mathrm{im}}=512$ in $\approx 2$ seconds. GPU-acceleration allows for simulating both more complex (realistic) physical models and larger image sizes (less stochasticity in likelihoods, see Section 2.5) than would be feasible with CPUs alone.

2.5 Likelihoods and Posterior Sampling

We evaluate the likelihood of a pCMD given a model with a binned Hess diagram and Gaussian statistics. With this approach, we create a binned 2D histogram of pixels in the pCMD and compare the relative number of counts in the two distributions, normalized to the total number of pixels. We choose bins of width $0.05$ magnitudes in each dimension, a width chosen to be roughly equivalent to the observational uncertainties, but we show in Section 3.2 that the resulting posteriors are fairly insensitive to this choice. We compute the log-likelihood $\mathcal{L}$ using the counts of data pixels $d_{i}$ and model pixels $m_{i}$ in each Hess bin, as:

[TABLE]

where the uncertainty $\sigma_{i}$ is approximated by the square-root of the number of Hess bin counts, added in quadrature:

[TABLE]

We apply a floor to the uncertainty to down-weight very rare bins in the Hess diagram, where differences in simulated and data image sizes may unintentionally bias the likelihood. This is one aspect of a general problem of evaluating likelihoods in Hess diagram space: how to handle data points when there are no (or very few) model points. For instance, a Poisson likelihood model would return zero likelihood if the model predicts zero pixels in a bin with even a single data pixels.

We add an additional likelihood term, corresponding to the difference in mean color $C$ and magnitude $M$ between the two distributions, with error of $0.05$ magnitudes:

[TABLE]

Without this term, two pCMDs which do not overlap at all ($d_{i}=0\;\forall\,i\;\mathrm{s.t.}\,m_{i}\neq 0$) are equally poor fits regardless of whether they are offset by an average of 1 or 10 magnitudes. The addition of this term gives slight preference to models where the center of the distributions roughly align with the data, without substantially affecting best-fit estimates since the relative magnitude of the term is quite small.

We sample the posterior using a new python package for dynamic nested sampling, dynesty (Speagle in prep). Nested sampling (e.g., Skilling, 2004; Feroz et al., 2009, 2013; Handley et al., 2015) is an approach similar to the commonly-used Markov Chain Monte Carlo (MCMC) technique, representing the posterior through a collection of samples from the distribution. Unlike MCMC, nested sampling efficiently computes the Bayesian evidence (also called marginal likelihood), allowing for principled model comparisons. Nested sampling algorithms also allow for sophisticated handling of multi-modal distributions. See the references above for more details on nested sampling.

Throughout this work, parameter estimates are reported as the median of the marginalized posterior probability function, and error-bars are reported as the $68\%$, equal-tailed credible interval, unless otherwise stated.

The pCMD likelihoods are stochastic, meaning that recalculating the likelihood multiple times for the same input parameters will result in a different log-likelihood. This presents a statistical challenge for any posterior sampling algorithm, and without proper accounting leads to an underestimate in uncertainties (see, however, the pseudo-marginal MCMC approach Andrieu & Roberts, 2009). We discuss this further in Appendix A, and detail a method for post-processing the results to recover reasonable posteriors.

Fitting a pCMD with models of size $\mathrm{N}_{\mathrm{im}}=512$ takes $\sim$100 GPU-hours for most cases. This is a speed-up of a factor of $7\times$ compared to the original code of Conroy & van Dokkum (2016), which required around 700 CPU-hours to fit a posterior using $\mathrm{N}_{\mathrm{im}}=256$. This was largely due to the need to combine the posteriors of 10 independent MCMC runs to overcome the stochasticity of the likelihoods. The larger simulated image sizes allowed by GPU-acceleration decrease the inherent stochasticity of the likelihoods, making combining multiple fits unnecessary.

3 A Suite of Mock Tests

3.1 Mock pCMD Models

To evaluate the capability of the code to recover model parameters, we run a series of mock tests where we fit models to simulated pCMDs generated from the code. In most cases, we fit models with the same physical components (i.e. same metallicity, SFH, and dust model) and hyper-parameters as used to generate the data. We model observations in the F814W ($\mathrm{I}_{814}$) and F475W ($\mathrm{g}_{475}$) bands of ACS, use a mock image size of $\mathrm{N}_{\mathrm{im}}=256$ and model image size of $\mathrm{N}_{\mathrm{im}}=512$, and assume no sky noise. Exposure times for F814W and F475W are set to 3235 and 3620 seconds, respectively, corresponding to the typical values of data from the PHAT survey (Dalcanton et al., 2012).

The fiducial model we frequently study, which we denote the "fiducial-$\tau$" model, has a fixed-width Gaussian MDF, tau-SFH, fixed-width log-normal dust screen, and the distance is allowed to vary (components M2+S3+E2+D2, see Section 2.2 for details). Unless otherwise specified, the free parameters used to generate the mock pCMDs are set to $[\mathrm{Fe}/\mathrm{H}]=-0.25$, $\log\mathrm{E}(\mathrm{B}-\mathrm{V})=-0.5$, $\mathrm{N_{pix}}=10^{2}$, $\tau=3$ Gyr, and $\mu_{d}=26.0$. The width of the MDF is set to $\sigma_{[\mathrm{Fe}/\mathrm{H}]}=0.2$, and the width of the Log-Normal dust model is set to $\sigma_{\mathrm{dust}}=0.1$.

We also study a "fiducial-nonparametric" model, where the SFH is fit with a 5-bin, non-parametric model (S5). The five SFH bins correspond to the following ages:

• •

SFH0: 1 Myr $-$ 100 Myr

• •

SFH1: 100 Myr $-$ 1 Gyr

• •

SFH2: 1 Gyr $-$ 3 Gyr

• •

SFH3: 3 Gyr $-$ 10 Gyr

• •

SFH4: 10 Gyr $-$ 14 Gyr

Unless otherwise specified, the same free parameters are used to generate mock pCMDs as above, with the exception of the SFH, which is a constant star-formation rate model with $\mathrm{N_{pix}}=10^{2}$.

In both models, we assume flat priors over all parameters. In most cases, we assume $[\mathrm{Fe}/\mathrm{H}]\in[-0.5,+0.25]$, $\log\mathrm{E}(\mathrm{B}-\mathrm{V})\in[-1,0]$, $\log\mathrm{N_{pix}}\in[2,5]$, $\tau\in[0.1,8.0]$ Gyr, and $\mu_{d}\in[22,26]$. In the case of the non-parametric SFH, we assume flat priors in $\log\mathrm{SFH}$, of width $\pm 1$ dex around the true underlying SFH.

3.2 Recovery of Non-Parametric SFHs

Figure 4 shows the posterior probability distribution for a fiducial-nonparametric model fit, and demonstrates that PCMDPy can simultaneously recover the input model parameters well. The metallicity, dust content, distance, and total $\mathrm{N_{pix}}$ (computed as the sum of all SFH bins) are all recovered to within $1\sigma$, as shown in the marginalized histograms. There is a notable degeneracy between $\mathrm{N_{pix}}$ and $\mu_{d}$, but the distance modulus is still constrained to within 0.1 dex.

The oldest bins of SFH (bins 3 and 4) are somewhat degenerate, but the total star formation older than 3 Gyr is fairly well constrained, as shown in the marginalized contours of those two bins. Figure 4 also shows the derived constraints on the SFR in each pixel, as a function of age. Compared to the original prior bounds allowed ($\pm 1$ dex in each bin), the model has recovered reasonably tight constraints on the SFH.

We also examine the recovered SFH for various input SFHs. The results are shown in Figure 5, for constant-SFR and $\tau=3$ Gyr models and increasing $\mathrm{N_{pix}}$ from $10^{2}$ to $10^{5}$. Each model is able to recover the underlying SFH, with the true SFH contained within the $68\%$ credible interval in nearly all cases. In constant-SFR models, the oldest ages of star formation have the largest uncertainties, while the relatively higher SFR in the $\tau$ models is easier to constrain precisely.

3.3 Evaluation of Model Choices

We use the suite of mock tests to evaluate the effect of various hyper-parameter choices and model families on the recovered parameters. For simplicity, these tests are performed using the fiducial-$\tau$ model, and the results are shown in Figure 6. These findings may not generalize fully to all models or regions of parameter space (especially higher $\mathrm{N_{pix}}$), and we caution all users of PCMDPy to think carefully about all hyper-parameter choices before fitting.

The first column of Figure 6 shows the effect of fitting distance as a free parameter (D1 vs. D2). We find we can recover the distance modulus to $\pm 0.1$ magnitudes. With the exception of $\mathrm{N_{pix}}$ (which is slightly degenerate with $\mu_{d}$), estimates of other parameters are just as well constrained regardless of whether distance is included in the fit.

The second column shows the effect of fitting for the width of the MDF and Dust distributions (Simple: M2+E2. Complex: M3+E3). We find the model is unable to constrain these widths, which leads to a small bias and inflated uncertainties in the estimates of other parameters. We therefore recommend against fitting for these width parameters, unless available information suggests a fairly informative prior is warranted.

In the third column, we vary the model image size $\mathrm{N}_{\mathrm{im}}$. Smaller image sizes ($\mathrm{N}_{\mathrm{im}}=256$) lead to inflated uncertainties, despite taking just as long to fit as the default size of $\mathrm{N}_{\mathrm{im}}=512$. The additional stochasticity in the likelihoods at small $\mathrm{N}_{\mathrm{im}}$ leads to more severe drops in sampling efficiency (see Appendix A). Larger image sizes ($\mathrm{N}_{\mathrm{im}}=1024$) take significantly longer to fit, requiring fewer nested sampling live points to converge in a reasonable time (fewer than 300-GPU hours), which results in under-estimated errors.

The final column shows various choices for the width of the Hess bins used to evaluate likelihoods. The bin size has no significant effect on the recovery of the posterior, although we note that this may not hold true at higher $\mathrm{N_{pix}}$, when the pCMD distribution is more compact.

In Figure 7, we study the ability of PCMDPy to model galaxies as a function of distance. At larger distances, the magnitude of surface-brightness fluctuations decreases (since $\mathrm{N_{pix}}$ is larger), making the models less sensitive to the underlying stellar populations. In addition, the same physical size region will correspond to far fewer pixels, and gradients in surface-brightness and stellar populations will make the modeling assumptions invalid over too many pixels.

To simulate these effects, we begin with a fiducial-$\tau$ model with $\mathrm{N_{pix}}=10^{2}$, modeled at distance of 1 Mpc, roughly equivalent to the disk of a galaxy like M31 (Conroy & van Dokkum, 2016, see also Cook et al. In Prep). The mock pCMD for this system is generated with $\mathrm{N}_{\mathrm{im}}=2048$, corresponding to a region 500pc on a side with total mass $\mathrm{M}_{\star}\approx 10^{8}\,\mathrm{M}_{\odot}$. We then approximate the effects of observing the same physical system at larger distance and correspondingly-higher $\mathrm{N_{pix}}$, while decreasing $\mathrm{N}_{\mathrm{im}}$ of the mock pCMD to keep the physical size and total mass constant. Five such mock datasets (from 1 to 100 Mpc) are fit to the same model (with simulated $\mathrm{N}_{\mathrm{im}}=512$).

The fits from this experiment are in excellent agreement with the true models until $\mathrm{D}\gtrsim 10$ Mpc ($\mathrm{N_{pix}}\gtrsim 10^{4}$), at which point the uncertainty rises sharply. But even out to 100 Mpc, the true parameters remain within the $68\%$ credible interval.

The pCMD method should therefore provide interesting constraints of stellar populations and distance out to large distances, and may serve as a useful complement to the existing SED-modeling technique.

To test the stability of the results, we run eight fits on different realizations of pCMDs generated with the same underlying model. The scatter in the estimates of each fit should fall within the typical uncertainty, or else the uncertainties are likely underestimated. Figure 8 shows the results of this test. For all parameters, the scatter (standard-deviation) in median estimates is significantly less than the typical uncertainty. This suggests that the model uncertainties are likely over-estimated by a factor of roughly $2\times$. It could be possible to correct for this with a more careful selection of $\mathcal{L}_{\mathrm{max}}$, defined in Appendix A, but we are comfortable with over-estimating our uncertainties, given the relatively loose method for dealing with the stochastic likelihoods detailed in the Appendix.

3.4 Model Mismatch Tests

We investigate the effect of fitting a pCMD with a different model than used to generate the data in Figure 9. First, we show fits where the physical model is incorrect. Incorrectly modeling the metallicity distribution function may lead to a minor systematic bias in $[\mathrm{Fe}/\mathrm{H}]$ and $\tau$, but for $\sigma_{[\mathrm{Fe}/\mathrm{H}]}{}$ of order a few-tenths of a dex, the effects are small and within the uncertainties.

Slight biases also arise from incorrectly modeling the SFH or the structure of the dust. We fit a model with $\mathrm{F}_{\mathrm{dust}}=1.0$ to data generated with $\mathrm{F}_{\mathrm{dust}}=0.5$, effectively ignoring that half of the stars should lie above the screen of dust. This results in an underestimate of the dust content, as well as bias in $\tau$ and $\mathrm{N_{pix}}$. Accurately accounting for the geometries of stars relative to the dust is therefore important to recovering accurate measurements.

We also show the effect of modeling a $\tau$-SFH with a non-parametric model. For comparison to the $\tau$ models, we compute the average age of the inferred SFHs, and convert to an effective $\tau$. The non-parametric model appears slightly biased towards preferring old ages (lower $\tau$). This results in slight positive bias in $\mathrm{N_{pix}}$ and $\mu_{d}$ in order to result in the same total luminosity. The age bias could be mitigated by increasing the number of SFH bins modeled, but this becomes computationally demanding. Incorrectly specifying the physical model can lead to subtle systematic biases, and it is not always easy to diagnose such a mismatch. One possible approach could be to compare the Bayesian evidence (computed through the nested sampling algorithm) for multiple model choices, or alternatively by studying the residual patterns in Hess diagram space. In future work, we additionally intend to study the effect of different SFH priors, such as those suggested by Leja et al. (2018).

Other possible mismatches between data and model include errors in the calibration of the observational data. If the exposure time of images is overestimated (the model shown overestimates the true exposure by $2\times$), the model greatly overpredicts the age parameter $\tau$, and most other parameters are also biased from their true values. When extracting pCMDs from observations, it is important to properly understand the observing conditions and whether multiple exposures were co-added to produce the final photometry.

The final two examples showcase the complex effects of poor PSF calibrations. If both PSF models are too wide (FWHM overestimated by $10\%$), best fit models are strongly biased in all parameters. Yet if only one PSF model is miscalibrated (here, the F814W filter), many of the biases go in the opposite direction. We investigate this effect more fully in upcoming work (Cook et al., In Prep) when we discuss the important practical challenges in selecting observational data to analyze with the pCMD technique.

4 Conclusions

We have presented a study of the application of pixel color-magnitude diagrams (pCMDs) for studying galaxies in the semi-resolved regime, as first introduced by Conroy & van Dokkum (2016). We developed the open-source, GPU-accelerated Python package PCMDPy, detailed its primary modeling assumptions and physical models implemented. We also presented the results of a suite of mock tests that demonstrate the potential for using PCMDPy to constrain important physical parameters of galaxies.

The main results of this work are as follows:

We have developed the new package PCMDPy for modeling pCMDs and inferring the physical properties of galaxies. It contains a number of improved physical models compared to the original code of Conroy & van Dokkum (2016), including a Gaussian metallicity-distribution function, a Log-Normal dust extinction model, and allows for simultaneous fitting of distance modulus. The code uses GPU-acceleration to allow the modeling of pCMDs with these complex models and larger simulated image sizes in significantly less time. Posteriors are sampled using the dynamic nested sampling code dynesty, and fitting a model requires $\sim 1/7^{\mathrm{th}}$ as much computational time as the original code. 2. 2.

We demonstrated that pCMDs are sensitive to the distance to a galaxy, and that our code can simultaneously recover accurate estimates of distance, SFH, metallicity, and dust content. 3. 3.

We show that highly precise measurements of distance and stellar populations should be recoverable for massive galaxies at least to 10 Mpc with HST-like resolution, and reasonable constraints are possible out to 100 Mpc. 4. 4.

The model is relatively robust against small biases in the assumed metallicity distribution function. Mischaracterized dust or SFH models can lead to small systematic bias in recovered parameters. Additionally, it is crucial to have good estimates of optical properties such as exposure time and well-calibrated PSF models. Future work will discuss other important considerations for selecting observational data to fit with PCMDPy. 5. 5.

The effects of the stochastic likelihood model present challenges for traditional sampling methods. The larger images that can be simulated using GPUs reduce this burden somewhat, but additional steps must be taken to account for the bias towards positive-likelihood fluctuations. Our chosen approach may result in overestimating the derived uncertainties by a factor of order $2\times$.

We thank Aaron Dotter and Jieun Choi for helpful discussions around the details of the MIST models. B.C. acknowledges support from the NSF Graduate Research Fellowship Program under grant DGE-1144152. This work is supported in part by HST-AR-14557. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

Appendix A Stochastic Likelihoods

The forward-modelling procedure for generating pCMDs is stochastic, and therefore so are our likelihoods. The simulated pCMD is a random realization of the underlying model, and therefore evaluating the likelihood of the same point in parameter-space multiple times will produce different results. Our tests indicate the distribution of likelihoods, given fixed model parameters, has a well-defined center but a long positive tail, essentially resulting from "overfitting" to the noise in the data. Larger simulated image sizes $\mathrm{N}_{\mathrm{im}}$ produce pCMDs which average over the many rare fluctuations, and therefore decrease the inherent stochasticity.

Any level of stochasticity poses significant problems for nested sampling or MCMC, as a positive fluctuation in likelihood will result in a model point being given larger weight than it truly deserves. In either MCMC or nested sampling, this has the detrimental effect of significantly decreasing sampling efficiency. The sampling algorithms will take longer to find a comparably good fit, resulting in decreased sampling efficiency and a cascading effect where most subsequent sampled points are also drawn from the positive-fluctuation tails.

Our tests have identified two other detrimental effects of nested sampling in a stochastic likelihood model. The weights assigned to the sampled points are biased towards only weighting the largest fluctuations, often resulting in only a single point being assigned a non-trivial weight. Since nested sampling integrates towards increasing likelihood, the final many points returned by the algorithm, and those assigned highest weight, almost all represent positive likelihood fluctuations, rather than actually better models.

Furthermore, the auto-stopping criteria used by many algorithms, including dynesty, relies on $\Delta\ln Z$, the estimated evidence remaining (e.g., Feroz et al., 2009). This value is biased high by the positive fluctuations, leading the algorithm to continue sampling longer than would be desired. Combined with the decrease in sampling efficiency, our dynesty fits often become stuck, taking thousands of likelihood calls to return a new sample888An example of this principle, for a simplified case of a stochastic likelihood, is included as a demonstration in the Dynesty respository:https://github.com/joshspeagle/dynesty/blob/master/demos/Examples--NoisyLikelihoods.ipynb. In practice, we stop the model fit after a fixed number of iterations. As discussed below, we find that we can later account for the bias in likelihoods, and when we do so the fits are almost always well converged by the $\Delta\ln Z$ criterion.

We show the effect of this likelihood bias in Figure 10. The blue line shows the measured log-likelihood of the dynesty samples from a mock test. The upturn at the final points shows extreme positive fluctuations: using the results straight from dynesty, the final point has $99\%$ of the weight, producing an implausibly peaked posterior. To demonstrate these points all represent positive likelihood fluctuations, we recomputed the likelihood of several of the points in parameter space, and show the measured distribution of likelihoods as error-bars.

It is clear that the likelihoods returned by dynesty are biased high from the average likelihood for each model. We can estimate what the maximum realistic likelihood should be, by computing many realizations of the "true" model used to generate the pCMD, and evaluating the fits. This is shown in the grey band. As should be expected, all of the mean likelihoods for the resampled points lie roughly within this band. We therefore have a means for estimating which likelihoods are trustworthy, and which are likely biased high due to fluctuations.

We use this insight to adopt a post-processing procedure to account for this likelihood bias. At the completion of each dynesty run, we recompute the likelihood of the best fit model 100 times, and then adopt the median value as a likelihood threshold, $\mathcal{L}_{\mathrm{max}}$. Any likelihoods above that threshold are capped to that value, and we recompute the weights and convergence statistics using the updated likelihoods.

We find that this approach is successful at eliminating the bias from likelihood fluctuations, by down-weighting the final handful of samples that are most subject to this stochastic bias. Figure 11 shows the mean and uncertainties in model parameters recovered as a function of this $\mathcal{L}_{\mathrm{max}}$. When $\mathcal{L}_{\mathrm{max}}$ is too high, the mean estimates vary sharply, as only one or a few sampled points are given all the posterior weight. Using our adopted $\mathcal{L}_{\mathrm{max}}$ ceiling, we are averaging over sufficiently many sampled points to recover stable means and reasonable uncertainties. When the remaining-evidence criterion is recomputed using the adjusted likelihoods, the fits are almost always extremely well converged by standard stopping criteria ($\Delta\ln Z\lesssim 0.01$).

This method for handling the stochastic likelihoods is a rough, heuristic approach, so the exact statistical uncertainties must be treated cautiously. The mock tests shown in this work are therefore important for validating that the approach recovers reasonable errors. As shown in Figure 8, the scatter in mean estimates between multiple independent trials is smaller than the typical uncertainty, indicating that if anything this relatively ad-hoc procedure is somewhat overly conservative in the model uncertainties.

Appendix B Cloud Computing and GPUs

The PCMDPy GPU-acceleration is written in CUDA, a proprietary language that only operates on Nvidia-brand GPUs. This limits the systems that can utilize the GPU-accelerated version of PCMDPy. Most personal computers and laptops (especially Macs) do not have Nvidia cards; Intel and AMD are common competitors, but CUDA code cannot be compiled to run on their architecture. However, the GPUs installed in most compute clusters hosted at research institutions are Nvidia, and cloud computing offers an alternative option for obtaining GPU-accelerated computational resources.

To make PCMDPy’s GPU capabilities more broadly applicable, we have also developed a Docker container architecture for fitting pCMDs in the cloud using AWS-Batch. Using this framework, pCMD models can be fit by renting GPU-enabled instances on-demand, with the benefit of having practically limitless scaling potential (large numbers of jobs can be run simultaneously). As of this writing, such instances are available for around \1$/GPU-hour (with faster GPUs available for higher cost, in nearly linear proportion to their speed-up). As PCMDPy typically takes around 100 GPU-hours to fit a model, this translates to around$$100$ per model fit. Full details and examples will be detailed later with the public PCMDPy code release paper. For a detailed example on using AWS in astronomical research, see Williams et al. (2018).