Measuring Star Formation Histories, Distances, and Metallicities with Pixel Color-Magnitude Diagrams II: Applications to Nearby Elliptical Galaxies
DDraft version April 23, 2020
Preprint typeset using L A TEX style AASTeX6 v. 1.0
MEASURING STAR FORMATION HISTORIES, DISTANCES, AND METALLICITIES WITH PIXELCOLOR-MAGNITUDE DIAGRAMS II: APPLICATIONS TO NEARBY ELLIPTICAL GALAXIES
B. A. Cook , Charlie Conroy , and Pieter van Dokkum Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA Astronomy Department, Yale University, 52 Hillhouse Ave, New Haven, CT 06511, USA
AbstractWe present spatially-resolved measurements of star formation histories (SFHs), metallicities, anddistances in three nearby elliptical galaxies and the bulge of M31 derived using the pixel color-magnitude diagram (pCMD) technique. We compute pCMDs from archival
HST photometry of M87,M49, NGC 3377 and M31, and fit the data using the new code
PCMDPy . We measure distances toeach system that are accurate to ∼ . The recovered non-parametric SFHs place reasonable ( ± dex) constraints on the recent (< 2 Gyr) star formation in M31 and NGC 3377, both of which showevidence of inside-out growth. The SFHs in M87 and M49 are constrained only at the oldest ages.The pCMD technique is a promising new avenue for studying the evolutionary history of the nearbyuniverse, and is highly complementary to existing stellar population modeling techniques. INTRODUCTIONThe evolution of galaxies, in particular the build-upof their stellar mass, is shaped by a combination ofmany physical processes including mergers, feedbackfrom supernovae and active galactic nucleii, and theaccretion of gas from the circumgalactic and intergalacticmedia. Constraining the impacts of these processes oftenrequires comparisons between observed galaxies and theoutputs of hydrodynamical simulations, which can studythese physical mechanisms in great detail (e.g., Hopkinset al. 2014; Vogelsberger et al. 2014). A galaxy’s starformation history (SFH) and chemical enrichment areparticularly sensitive to the relative contributions of theseprocesses.One common method for measuring stellar popula-tions and SFHs in galaxies is resolved-star photometry(e.g., Dolphin 2002; Weisz et al. 2011; Lewis et al. 2015;Williams et al. 2015), which compares the colors andmagnitudes of individual stars to stellar evolution models.Due to the limited angular resolution of even the bestoptical telescopes, SFHs measured with the resolved-starphotometry technique are primarily limited to galax-ies in the Local Group and slightly more distant dwarfgalaxies, due to crowding. Because young, massive main-sequence stars can be identified even against a crowdedbackground, constraining more recent epochs of starformation with resolved stars is notably easier than mea-suring ancient SFHs, which requires resolving stars as [email protected] faint as the oldest main-sequence turnoff for the resultsto be robust against systematic uncertainties (Schulzet al. 2002; Williams et al. 2017).An alternative method, spectral energy distribution(SED) modeling, compares broadband photometry orspectra of the integrated light of (typically) the entiregalaxy against stellar population synthesis models torecover stellar populations and SFHs (e.g., Walcher et al.2011; Conroy et al. 2013). The advent of integral-fieldunits (IFUs) allows for spatially-resolved measurementsof nearby galaxies with SED-modeling techniques (e.g.,Smith & Hayward 2018). But SFHs derived from SED-modeling are most sensitive to the levels of ongoing starformation, as the light from old, low mass main-sequencestars is usually overwhelmed by the light of a few rare,evolved, young stars (Papovich et al. 2001; Marastonet al. 2010; Pforr et al. 2012; Sorba & Sawicki 2015; Lejaet al. 2018). The high signal-to-noise required to preciselymeasure stellar populations using IFUs also often limitsspatially-resolved measurements to the brightest innerregions of galaxies.These two techniques can be categorized within a com-mon framework by considering the typical number ofstars per resolution element, a quantity denoted N pix (vanDokkum & Conroy 2014; Conroy & van Dokkum 2016).The unresolved measurements used in SED-modelingtypically have N pix (cid:38) O (10 ) , while robustly measuringSFHs by fully resolving the oldest main-sequence turnoffrequires N pix (cid:46) O (10 − ) .In between these two observational regimes lies theso-called semi-resolved regime , defined very roughly as a r X i v : . [ a s t r o - ph . GA ] A p r O (10 ) (cid:46) N pix (cid:46) O (10 ) . All but the brightest individ-ual stars cannot be detected in semi-resolved photom-etry, but surface-brightness fluctuations from Poissonsampling of rare, bright stars in each pixel are significant.Massive galaxies viewed with the Hubble Space Telescope( HST ) are typically semi-resolved from ∼ Mpc out tonearly 100 Mpc.A method for measuring stellar populations and SFHsin this regime, known as the pixel color-magnitude dia-gram (pCMD) technique, was first introduced in Conroy& van Dokkum (2016) as a way to model semi-resolvedphotometry. The authors demonstrated that the proper-ties of pixel-to-pixel surface-brightness fluctuations aresensitive to the underlying stellar populations and SFH,and made a first application of the method to the bulge ofM31. Several other works have also studied stellar pop-ulations through the distribution of pixels in CMD space,including either qualitative comparisons of pCMDs be-tween galaxies (Bothun 1986; Lanyon-Foster et al. 2007;Lee et al. 2017) or quantifying the dispersion of thepCMD (Lee et al. 2018). However, these works focus onthe internal variation of stellar populations across pixels,rather than using the surface-brightness fluctuations asa probe of the stellar populations themselves.Cook et al. (2019) presented
PCMDPy , a new Pythonpackage developed to implement the pCMD technique,with additional levels of model complexity and the abil-ity to simultaneously measure distances. Mock testsof
PCMDPy demonstrated the potential for the pCMDtechnique to measure spatially-resolved star formationhistories (SFHs) and abundances in systems to distancesof at least 10 Mpc, and possibly as far as 100 Mpc.The pCMD technique was also shown to be capableof simultaneously measuring distances to galaxies. Asdescribed in Cook et al. (2019), the pCMD technique canbe considered a generalization of the surface-brightnessfluctuations (SBF) distance technique (Tonry & Schnei-der 1988). The SBF metric has been applied to ellipticalgalaxies and large bulges of spiral galaxies (Ferrareseet al. 2000), but absolute calibration is fairly dependenton uncertain, evolved phases of stellar evolution (Tonryet al. 2001), and the effects of varying stellar populationscan be significant.We present here the first application of the pCMD tech-nique to galaxies outside the Local Group, and the firstmeasurement of distances from pCMDs. In Section 2,we describe a procedure for cleaning, reducing, and align-ing archival
HST images in order to compute pCMDsthat can be analyzed with
PCMDPy . In Section 3, we out-line the models used to fit these data. The results arepresented in Section 4, including the measured star for-mation histories (4.2), distances (4.3), metal abundanceand dust extinction (4.4). We discuss future considera-tions for the pCMD technique in Section 5, and conclude in Section 6. HST
DATA REDUCTION2.1.
Data Collection
We compute pCMDs using archival photometry fromthe Hubble Space Telescope’s Advanced Camera for Sur-veys (
HST -ACS). We choose three elliptical galaxieswith multiple exposures in at least two wide-band filters,and total exposure time in each filter of at least 700seconds. The galaxies studied are M87, M49, and NGC3377. We also study five regions in the bulge of M31, us-ing pre-drizzled data from Brick 01 of the PHAT survey(Dalcanton et al. 2012). We adopt distance modulus ( µ d )values from the literature, as shown in Table 1, which areused to compute radial distances within a galaxy, and forcomparison to our measured distances. The archival datacollected for each galaxy are summarized in Table 2.In most galaxies, we use photometry in only two filters:a red ( I or z ) and a blue ( g ) filter. The exceptionis M87, where V photometry with significantly longerexposure time than the g data is also available. TheM87 results use the I − g color, but we find identicalresults, within the uncertainties, when using the V photometry in place of g .2.2. Combining and Aligning Exposures
We align each set of individual .flc exposures using
TweakReg from the
AstroDrizzle package, and thencombine exposures with the drizzling tool (Fruchter et al.2010; Gonzaga et al. 2012). This increases the effec-tive exposure time of the photometry, removes cosmicrays, and corrects for non-linear optical distortions. The
TinyTim package (Krist et al. 2011) does not providepoint-spread function (PSF) models that account fordrizzling, so we use the lanczos3 kernel, which mostfaithfully preserves the properties of the original PSF(Jee et al. 2007). We disable sky subtraction, choosinginstead to forward model the addition of sky backgroundduring the simulation step (see Section 2.5). Regardless,several of our sources are sufficiently extended that thedefault
AstroDrizzle sky subtraction algorithm wouldsignificantly overestimate the background. The defaultsare used for all other
AstroDrizzle parameters.Modeling of pCMDs assumes that the images in eachfilter are extremely well aligned. Yet after the initialdrizzling procedure, we find that offsets of as much asa few pixels can remain between filters. We thereforerealign and re-drizzle the final images to a commonreference frame, following Section 7.5 of Gonzaga et al.(2012). After this second alignment, the images of eachgalaxy in separate filters are found to be aligned to betterthan 1/10 th of a pixel, according to TweakReg .2.3.
Data Cleaning and Source Removal
Galaxy µ d (mag) D (Mpc) Method Dist. Reference [Z / H] [Z / H] ReferenceM31 24.44 0.77 TRGB Conn et al. (2016) 0.23 Conroy & van Dokkum (2012)NGC 3377 30.2 10.9 TRGB Lee & Jang (2016a) 0.13 Conroy & van Dokkum (2012)M87 30.9 15.1 TRGB Lee & Jang (2016b) 0.17 Conroy & van Dokkum (2012)M49 31.1 16.8 SBF Blakeslee et al. (2009) 0.30 van Dokkum & Conroy (2014)
Table 1 . Benchmark literature distances to the galaxies and measured metallicities. We adopt TRGB distances when available,then use distances derived from the SBF technique. We assume these are the true distances, enabling an assessment of thedistance measurements with the pCMD technique. Radial distances within a galaxy are determined assuming the distancesgiven here. The metallicities listed include variable abundance ratios and computed at various radii, but should be roughlycomparable to those measured in this work.Galaxy Filter Exposures Exp. Time (s) Proposal ID Sky (counts)M31 F814W (I) · · · · · ·
Table 2 . Summary of photometric data collected from
HST archives. We use drizzled mosaics from the PHAT survey for M31,and so do not reduce individual exposures. The final column lists the estimated sky background, computed using the ACSExposure Time Calculator, and is in units of counts per pixel, integrated over the entire exposure (see Section 2.5).
The mock tests in Cook et al. (2019) showed thatimproperly estimating the exposure time can significantlybias the resulting fits, because the effects of photon noiseon the pCMD distribution are challenging to disentanglefrom the true surface-brightness fluctuation signal . Wetherefore mask any pixels that have an effective exposuretime (as computed by AstroDrizzle ) lower than ofthe total in either band. This removes pixels that werenot contributed to by all individual exposures. Commoninstances include pixels affected by cosmic rays, hotpixels, or regions where the exposures did not overlap,such as across the ACS chip gap. This results in maskingapproximately of all pixels.We mask all detectable sources in the images (globu-lar clusters, foreground stars, and background galaxies)which will have significantly different properties in pCMDspace than the diffuse field star populations in the targetgalaxy. We first use Source Extractor (Bertin & Arnouts1996; Barbary & Kyle 2016) to subtract a smoothed To account for this, the SBF distance technique isolates thetrue signal by analyzing the fluctuations in Fourier space. model of the galaxy and then automatically identify andmask the majority of sources. We then mask by-eye anyresidual sources, such as obvious globular clusters neargalaxy centers, and stray light from extended backgroundgalaxies and bright foreground stars.2.4. pCMD Extraction
The final step in extracting pCMDs from data is toidentify regions where the standard assumption of model-ing pCMDs – that the typical number of stars in a pixel( N pix ) is roughly constant – is at least approximatelyvalid. To this end, we focus on elliptical galaxies (andthe bulge of M31) in this work, as they typically showfar less substructure than spiral galaxies. Yet all galax-ies exhibit significant radial surface-brightness gradients,making it impractical to apply the pCMD technique tothe entire galaxy.To overcome this limitation, we extract pCMDs fromthin, roughly elliptical regions in each galaxy, withinwhich the average surface-brightness does not vary sig-nificantly. This yields the added benefit of allowing usto measure spatially-resolved SFHs and metallicities asa function of radius. In selecting such regions, however,there is a balancing act required between selecting slicesthat are too wide (such that surface-brightness gradientsdistort the pCMD distribution) and those that are toothin (such that too few pixels remain to be analyzed).For this reason, we are unable to extract valid pCMDsin the central regions of each galaxy (roughly within10"), and instead focus on the outer regions, as describedbelow.We compute elliptical slices aligned with the galaxiesusing SAOImage DS9 (Joye & Mandel 2003), by creat-ing highly smoothed contours (smoothness ≈ ) andconverting them to polygonal regions. The contours arescaled logarithmically, and spaced such that the typicalsurface-brightness gradient across each region is signifi-cantly less than 0.1 magnitudes. In practice, this resultsin regions of roughly 20 pixels in width. In each galaxy,we select three regions (five in M31), equally spaced inradius, and further divide these into four quadrants tominimize the effects of any large-scale azimuthal asym-metry. We primarily focus our analysis on one quadrantfrom each radius, although we confirmed that the resultsbelow are consistent for each quadrant.Any remaining large-scale variations in surface-brightness (either due to azimuthal variations or overly-wide annuli) in the extracted regions will result in broad-ened pCMD distributions. Mock tests indicate that theprimary impact on the inferred parameters is an under-estimate of N pix (and hence distance) to match the scaleof surface-brightness fluctuations, although the magni-tude of the bias is challenging to estimate and dependssignificantly on N pix .We compute pCMDs by first converting the instru-mental fluxes to apparent magnitudes, using zero-pointscomputed using the pysynphot package (Lim et al. 2015).We extract separate pCMDs from each of the regions,and remove the masked pixels as described above. Ex-ample pCMDs from three regions of M87 are shown inFigure 1. 2.5. Background Estimation
Background (sky) noise has a significant impact on theobserved pCMDs, and must be considered carefully. Ifwe consider sky noise to be a constant background in eachfilter, then its primary effects are to 1) increase the aver-age flux in each pixel, and 2) change the level of photonnoise measured at the CCD. Accurately modeling photonnoise is crucial to modeling pCMDs, because it can easilybe mistaken for intrinsic surface-brightness fluctuations(Cook et al. 2019). Therefore, we add a synthetic skysignal (a constant count-per-pixel rate in each filter) intoour simulated pCMDs, rather than subtracting the skybackground from our photometry.Several of our sources, especially M87 and M49, arewidely extended on the sky, such that the entire
HST field g - I I M87
Figure 1 . Example pCMDs extracted from quadrants of theM87 photometry, at the three radii given. Contours showthe bounds within which 39%, 87%, 99%, and 99.9% of thepixels lie (the σ , σ , σ , and σ contours, respectively). Thedecrease in average flux and increase in variance with radiusdemonstrates the effect of decreasing N pix . of view is contaminated by emission from their diffuseenvelopes. This makes automated methods for estimat-ing the sky, such as those provided in AstroDrizzle ,impractical, as the background would be overestimatedby including emission from the source.We instead estimate the sky background using
HST ’sonline Exposure Time Calculator (ETC) . Given theobserved filters, location of the source, and date of obser-vation, ETC computes the expected background countrate from zodiacal light, which we convert to counts-per-pixel integrated over the entire exposure. The skysignal we add into each simulated filter are listed inTable 2. These estimates also include the backgroundcontribution from Earthshine, using the models providedby the ETC and available pointing data from each ob-servation. In most cases the Earthshine background isvery small compared to the zodiacal light. We validateour background estimates in NGC 3377, which is far lessextended than M87 and M49, by computing the medianflux in the corners of the image farthest from the source.The values agree with our estimates to within . Thediffuse extragalactic background is expected to be lessthan a additional contribution (Zackrisson et al.2009), and is not included.Noise from dark current can also have a significanteffect on the level of photon noise, in some regions con-tributing as much as 10 % of the total flux. Dark currentis already subtracted from the .flc exposures, and wemask hot and warm pixels in the data cleaning stageabove. To account for the average dark current of the re-maining pixels, we add the ETC estimated dark current(0.0127 counts per pixel per second) to our simulated http://etc.stsci.edu/etc/input/acs/imaging/ images prior to applying Poisson photon noise, and thensubtract the same values from the results. METHODS3.1.
Overview of
PCMDPy
We analyze the data using
PCMDPy , a Python codedeveloped to fit pCMDs and infer the Bayesian poste-rior over physical parameters. We provide here a briefoverview of the code, while full details can be found inCook et al. (2019).Pixel Color-Magnitude Diagrams with Python(
PCMDPy ) is a GPU-accelerated python code for inferringphysical parameters of semi-resolved galaxies throughfitting observed pCMDs to synthetic models. It relies onthe key assumption of pCMD modeling, that the starsin each pixel are drawn via a Poisson process from thesame underlying population, with an average number ofstars per pixel given by a parameter N pix . The pCMDssimulated by PCMDPy therefore are only applicable tosmall regions of a galaxy where this assumption is valid.Simulating a pCMD with
PCMDPy begins with an as-sumed model for the metal abundances and SFH (withthe total stars formed summing to N pix ) of a region,and an initial mass function (by default, Salpeter 1955).Stars from this model are randomly populated into eachsimulated pixel of an image, and the flux in each pixelis computed using the MIST stellar evolution models(Choi et al. 2016). We apply dust attenuation, adjust forthe modeled distance to the source, and finally includeobservational noise such as sky background, PSF convolu-tion, and photon noise according to the properties of thetelescope and the simulated exposure time. The fluxesin each pixel are converted to magnitudes, resulting in apixel color-magnitude diagram.The primary free parameters of the physical modelspecify the metal abundance, star formation history, dustattenuation, and distance, with multiple options for eachmodel available in PCMDPy . The likelihood, or agreementbetween data and modelled pCMDs, is computed bybinning the pixels into a Hess diagram and comparing therelative counts. The posterior distribution is estimatedusing the nested sampling code dynesty (Speagle 2019).3.2.
PCMDPy models
Three model changes have been implemented relativeto that described in Cook et al. (2019). The first, asdescribed in Section 2.5, is the addition of dark currentnoise, which is added prior to simulating the Poissonphoton noise and then subtracted.Secondly, we have removed the sub-pixel PSF model,and instead apply a single PSF to the entire image,because subsequent tests indicate the particular sub-pixel PSF model described in Conroy & van Dokkum(2016) and Cook et al. (2019) is insufficient for modeling the desired effects. We refer the reader to the discussionin Section 5.3.Finally, we have added an additional metallicity distri-bution function (MDF) model that replicates the MDFof a "closed box" evolution model (Binney & Merrifield1998). The closed box MDF has a similar shape toa normal distribution, but benefits from having only asingle free parameter, as the width scales with the meanmetallicity. The closed box model may not, of course,be a perfect representation of the shape of the MDF ina narrow elliptical annulus of an ETG, and future workshould consider whether better physically-motivated dis-tributions, such as the "leaky-box" model, result in anysignificant changes to the derived physical properties.As our default model, we fit the pCMDs with a 5-binnon-parametric SFH, a single metallicity, and single dustextinction screen, and allow the distance modulus tovary (Models M1+S5+E1+D2, from Cook et al. (2019)Table 1), for a total of 8 free parameters. We assume flatpriors over all parameters in the model, including [Fe / H] (from -1.0 to 0.5), log E(B − V) (from -2.0 to -0.5), and µ d ( ± magnitudes around the literature distance). Foreach galaxy, we visually identify an approximate N pix through comparisons to simulated pCMDs with a 10 GyrSSP. We then assume a flat prior in the 5 star formationhistory parameters, with limits of ± dex around a τ = 3 Gyr model with that approximate N pix .Model images have N im = 512 ( pixels), and weassume a Salpeter IMF (Salpeter 1955). Posteriors arefit with the dynesty dynamic nested sampling algorithm(Speagle 2019), using 400 live points and uniform sam-pling within a multi-ellipsoidal boundary. Each fullmodel takes around 120 GPU-hours to fit, executed onNVidia Tesla K20xm chips.In addition to the default model described above, westudy two additional models. One assumes a fixed dis-tance, equal to the literature value given in Table 1, whilethe other replaces the single metallicity assumption witha closed box MDF.Throughout this work, parameter estimates and error-bars are reported as the median and credible intervalof the marginalized posterior distribution.3.3. Likelihood Model and Statistical Convergence
As discussed in Cook et al. (2019), we implement apost-processing correction to the final dynesty weights toaccount for biases arising from the stochastic likelihoodfunction. Given the very wide priors we assume here,we find the stochastic effects and decreased samplingefficiency to be much more significant than in earliermock tests. We found that in most cases, the fits did notstatistically converge. The sampling efficiency declinedso significantly that it took several thousand likelihoodcalls to produce each of the final few hundred sampledpoints.Even using the post-processing correction described inCook et al. (2019) often resulted in only a few represen-tative samples (with non-negligible posterior weights).We therefore applied a more liberal correction, lower-ing the log-likelihood ceiling until the nested samplingconvergence metric ( ∆ ln Z ) was less than . . Thisresults in weighting more evenly over many more of thefinal sampled points. We discuss this in greater detail inSection 5.1, but note that the errorbars presented heremay be overestimated. RESULTS4.1.
Overview
The final output of our modeling procedure are sam-ples representing the 8-dimensional posterior probabilitydistribution. A convenient way to visualize this distribu-tion is through a "corner plot", showing the marginalizedposterior distribution between all sets of two parameters.We show examples from Region E (1 kpc) in M31 and Re-gion C1 (6 kpc) in M87, in Figures 2 and 3, respectively.In each figure, we also show the inferred cumulative SFH,with the confidence region.In the outer M31 region, we see evidence for most ofthe same parameter degeneracies as in the mock tests ofCook et al. (2019). The distance is strongly degeneratewith N pix , and there is also a notable dust-metallicitydegeneracy. There are also strong correlations betweenthe star formation parameters, especially at the twooldest ages (SFH3 and SFH4), but together the overallamount of old star formation is well constrained.The M31 region shows strong evidence for relativelyelevated (more than a τ = 5 Gyr model) levels of starformation through the past ∼ Gyr, before apparentlyundergoing a quenching period. We see evidence of asharp decrease (by more than an order of magnitude)in the star formation rate between 2 and 0.3 Gyr ago.The old star formation is less well constrained in M87,but is marginally inconsistent with very rapid and earlyquenching (such as a τ = 1 Gyr model). We discussthis further in Section 4.2.In both cases, we are only able to place upper limits tothe amount of star formation in the youngest two bins.The dynesty fits do constrain the upper limits, but thetail of allowed star formation stretches to the lower limitsof our flat priors. Despite the large luminosities and veryblue colors of young, massive stars, they are so rare thatthey will populate only a very small number of pixels,and as such the data are not constraining between verylow levels and none at all. As such, even significantlyexpanding the prior volume will not result in a convergedlower limit. We show the upper limit in such cases,and discuss this further in Section 5.5. The models are able to replicate well the input datain each of the galaxies studied, and examples are demon-strated in the residual plots shown in Figure 4. ThepCMD distributions of data and best-fit models are over-laid, showing a high degree of similarity.4.2.
Spatially-Resolved Star Formation Histories
Among the primary results from the pCMD modelingare measurements of the spatially-resolved star formationhistory in each of the galaxies. The posterior distribu-tions of cumulative star formation, as a function of radius,are shown for the four galaxies in Figure 5.Conroy & van Dokkum (2016) previously applied thepCMD technique to derive the cumulative SFH in M31,and we show comparisons to their results in Figure 5. Atthe oldest ages, our measured SFHs agree well with theirresults, as we find evidence for increased levels of starformation in the outer regions, and evidence of earlierquenching in the inner ∼ . kpc. While the oldest binsof star formation agree within the uncertainties withConroy & van Dokkum (2016), we emphasize that webelieve the only an upper limit on the youngest periodsof star formation can be constrained.In NGC 3377, we see evidence for significant gradientsin the SFH. The outer region shows strong evidence ofelevated star formation, with as much as of thestars formed within around the last 300 Myr . Theinner region, by comparison, appears to have quenchedmuch earlier, with significantly less than of its starsforming in the last Gyr.In contrast to the relatively well constrained historiesin M31 and NGC 3377, the SFHs in M87 and M49 are,within the uncertainties, consistent at all radii. The old( (cid:38) Gyr) star formation in both systems is consistentwith τ models, with exponential timescales between − Gyr. Only upper limits can be placed on the youngeststar formation, but the limits are inconsistent with thelevels of recent star formation expected from a τ ≈ Gyrmodel. These findings, that old SFHs are consistent with τ ≈ Gyr while the recent SFHs are not, may indicatethat τ models are too inflexible to model realistic SFHs(e.g., Carnall et al. 2018; Leja et al. 2019).In Figure 6, we compare SFHs between galaxies, over-laying the posterior distributions at similar radii. Wedo not have overlap between all four galaxies at anyindividual radius.In the inner ∼ of its stars in the last ∼ Gyr. Yet at 4 kpc,NGC 3377 has a much younger SFH than M49 and M87.At all radii, the star formation histories of M49 and M87 although this is sensitive to SFH model, see 7 [Fe/H] = 0.24 +0.100.11 . . . . . l o g E ( B - V ) log E(B-V) = 1.33 +0.310.37 . . . . . l o gS F H log SFH0 = 1.72 +0.620.65 . . . . . l o gS F H log SFH1 = 0.64 +0.460.67 . . . . l o gS F H log SFH2 = 1.63 +0.220.19 . . . . . l o gS F H log SFH3 = 2.06 +0.270.88 . . . . . l o gS F H log SFH4 = 2.54 +0.160.26 . . . . . l o g N p i x log N pix = 2.70 +0.070.07 . . . . [Fe/H] . . . . . d . . . . . log E(B-V) . . . . . log SFH0 . . . . . log SFH1 . . . . log SFH2 . . . . . log SFH3 . . . . . log SFH4 .
40 2 .
55 2 .
70 2 .
85 3 . log N pix . . . . . d d = 24.67 +0.150.15 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r PosteriorPrior
Figure 2 . Recovered posterior probability distribution for M31 Region E, at a radius of 1.2 kpc. The σ , σ , and σ contoursare shown. The literature distance modulus is shown as a red line. The SFH parameters have units of stars-per-pixel, and N pix is computed from their sum. The upper-right panel shows the distribution of cumulative star formation history ( credible interval). The allowed prior region is shown as a dotted box, and dot-dashed lines show comparisons to τ -SFHs. Wecannot constrain lower limits to the star formation rates in the youngest two bins, but there is evidence of a decrease in the starformation rate of more than an order of magnitude over the log age interval from 9.25 to 8.5. [Fe/H] = 0.25 +0.380.33 . . . . . l o g E ( B - V ) log E(B-V) = 1.25 +0.410.49 . . . . . l o gS F H log SFH0 = 0.07 +0.850.83 . . . . . l o gS F H log SFH1 = 1.37 +1.091.07 . . . . . l o gS F H log SFH2 = 2.54 +1.081.39 . . . . . l o gS F H log SFH3 = 3.79 +0.471.21 . . . . . l o gS F H log SFH4 = 4.34 +0.250.98 . . . . . l o g N p i x log N pix = 4.51 +0.150.23 . . . . . [Fe/H] . . . . . d . . . . . log E(B-V) . . . . . log SFH0 . . . . . log SFH1 . . . . . log SFH2 . . . . . log SFH3 . . . . . log SFH4 . . . . . log N pix . . . . . d d = 31.35 +0.520.54 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r PosteriorPrior
Figure 3 . Same as Figure 2, for Region C1 of M87, at a radius of 6 kpc. g - I I M31 (1.2 kpc)DataModel g - I I = -1.06e+03 g - z z NGC 3377 (4 kpc) g - z z = -3.39e+02 g - I I M87 (6 kpc) g - I I = -8.04e+01 ± L o g - L i k e li h oo d p e r b i n Figure 4 . Comparisons between example datasets and their corresponding best-fit single metallicity models. The left columnshows the contours of the data (black) and model (green) pCMDs as in Figure 1. The right column shows the likelihood grid asa heatmap in Hess diagram space, multiplied by ± to show over- (under-)densities of model points in red (blue). The datapCMD is overplotted for reference. The total log-likelihood of the model is given in the bottom-right corner. The pCMDs andresiduals for M49 are qualitatively similar to M87, and so are not shown. The pCMD contours and residuals for the closed boxMDF models are also visually similar to the single metallicity models. are generally consistent, within the uncertainties.It is important to properly understand the effects ofdifferent model assumptions on the inferred stellar pop-ulations. We show, in Figure 7, the changes in derivedstar formation histories in M31 (Region E) and NGC3377 (Region C1) due to different assumed physical mod-els. These two regions were the only two with noticeablechanges in inferred SFH between models.When distance is held fixed, the inferred histories ofstar formation are remarkably consistent, especially giventhat the distance measured in the NGC 3377 regionshown is underestimated by around (see Section 4.3).The only notable offset is the second oldest SFH bin inNGC 3377, where the star-formation decreases slightly.Assuming a closed box MDF distribution rather thana single metallicity does slightly alter the inferred SFHs, suggesting somewhat earlier quenching in M31 and amarginally smoother SFH in NGC 3377.In M31 and NGC 3377, the spatially-resolved SFHsindicate a general pattern of more recent star formationwith increasing radius. This is further demonstrated inFigure 8, where we show the fraction of mass formedwithin the last 2 Gyr as a function of radius. Becauseof the small differences between the inferred SFHs bymetallicity model, we show the trends for both models,yet the results are similar in both cases.Despite the large uncertainties, the outer regions ofNGC 3377 and M31 both show evidence for more recentbuildup of mass: the outer M31 field formed more than of its stars in the last 2 Gyr, and as much as of the stars in the outer NGC 3377 region formed asrecently. These results are in general agreement with0 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r M31 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r NGC3377 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r M87 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r M49
Figure 5 . The cumulative SFH (with intervals) for three regions in the M31 bulge ( upper-left ), NGC 3377 ( upper-right ),M87 ( lower-left ), and M49 ( lower-right ). Prior boundary, upper limits, and example τ SFHs as in Figure 2. For clarity, we plotpoints with small offsets in age and replace shaded intervals with error bars. In the upper-left panel, we show comparisonswith Conroy & van Dokkum (2016), who first applied the pCMD technique to derive SFHs in M31, and our measurements agreewithin the uncertainties, although we find the measurements of young star formation to be only upper limits. Both M31 andNGC 3377 show evidence for longer periods of star formation in the outer regions, whereas the inner regions likely quenchedearlier. the inside-out formation scenario (White & Frenk 1991;Roškar et al. 2008; Rix & Bovy 2013; Patel et al. 2013).The young star formation in M87 and M49 is too poorlyconstrained to shown any indication of inside-out growth.4.3.
Distances with the pCMD Technique
In addition to measuring spatially-resolved SFHs, Cooket al. (2019) demonstrated the capability of measuringdistances with the pCMD technique. We recover a sepa-rate measurement of the distance modulus to a galaxyfrom each region. In Figure 9, we show the posteriordistribution of these distances, relative to the literaturevalues, for both the single metallicity and closed boxMDF models.In the single metallicity case, we find good agreement to within ± in all cases, with the exception of theouter region (C1) of NGC 3377 (see Figure 9), as thedistance estimate is biased low by around 0.6 magnitudes.The best-fit metallicity in this region (see Figure 11) isalso significantly higher than the other regions in NGC3377. We suspect this may be caused by a moderateunderestimate in the assumed sky background, whichthe models may compensate for by preferring a closerdistance to the source and which would affect the fainteroutskirts more significantly.In the case of the closed box abundance model, themeasured distances appear to be systematically biasedhigh by on average, but we do note the distances toNGC 3377 are in better agreement. The typical distances1 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r R 1.5 kpc
M31 BulgeNGC3377 log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r R 4 kpc
M87M49NGC3377
Figure 6 . Comparison of the cumulative SFHs of galaxies atsimilar physical radii.
Top:
Regions approximately 1.5 kpcfrom the center of M31 (1.2 kpc) and NGC 3377 (1.7 kpc).
Bottom:
Regions around 4 kpc from the centers of NGC3377(3.9 kpc), M87 (4.7 kpc) and M49 (4.4 kpc). to M87 and M49 are likewise biased high, but are consis-tent with the true values within the uncertainties. Weattribute the overestimated distances to the fact that theclosed box model is known to overestimate the width ofthe metallicity distribution functions in observed galaxies(e.g., Holmberg et al. 2007).We summarize the distance estimates to each of thefour galaxies in Figure 10, where we compare the dis-tances derived with the pCMD technique to literaturedistances measured with the TRGB and SBF techniques.We take the pCMD method estimate to be the weightedaverage of the individual measurements from each region.Both the TRGB distances and the SBF distances aresubject to systematic uncertainties of ∼ . magnitudes(not shown in Figure 10). In the case of TRGB, thisarises from the uncertainty in the TRGB magnitude (Lee& Jang 2016b), while for SBF it is due to calibration log age (yr) C u m u l a t i v e S F H = G y r = G y r = G y r M31, 1.2 kpc
Default ModelFixed DistanceClosed Box MDF log age (yr) C u m u l a t i v e S F H = G y r = G y r NGC3377, 3.9 kpc
Default ModelFixed DistanceClosed Box MDF
Figure 7 . The cumulative SFHs as a function of modelassumed. We show the outer regions studied in the M31bulge ( upper ) and NGC 3377 ( lower ), the only two datasetswith noticeable differences between the measured SFHs. relative to the Cepheid distance metric (Blakeslee et al.2009). The pCMD technique would likely be subject to asimilar level of systematic uncertainties as the TRGB, asit relies on well-calibrated luminosities of the brighteststars.The average distances measured with the pCMD tech-nique agree very well with the existing methods, althoughagain we note the potential systematic differences be-tween the two metallicity models. Including this possiblebias, the distances measured with pCMDs appear reliableto around − , and the pCMD technique is there-fore very useful for simultaneously making approximateestimates of distances and stellar populations of galaxieswithin 100 Mpc. Improving the precision of measureddistances will likely require significant improvements tothe likelihood model (see Section 5.1).4.4. Metal Enrichment and Dust Radius (kpc) C u m u l a t i v e S F H ( < G y r ) M87M49NGC3377M31
Figure 8 . Spatially-resolved young (< 2 Gyr) mass fraction. From each SFH distribution, we derive the fraction of mass formedwithin last 2 Gyr. The solid symbols show results from the single metallicity model, while open symbols show ages from theclosed box MDF models. Upper-limits are shown where the star-formation history over the past 2 Gyr was not constrained.M31 and NGC 3377 show hints of inside-out growth (more recent star formation at larger radius), but this is not constrained inM49 and M87.
The
PCMDPy models also derive estimates of the metalenrichment, in terms of the iron-abundance [Fe / H] .Figure 11 shows the metallicity derived at each radiusin the galaxies studied, although we emphasize that themetallicity is at best poorly constrained in most regions.This is primarily due to the strong degeneracy betweendust and metallicity (see Figures 2 and 3), which themodels are not able to easily distinguish between at thishigh N pix . Given the large uncertainties on measuredmetallicity, we can place no significant constraints onabundance gradients at the radii studied.The figure also includes the literature [Z / H] metallicities given in Table 1, which are notspatially resolved. The measured metallicities in M31agree very well with the literature. The metallicitiesmeasured in the other three galaxies are generally toolow, although they are within the large uncertainties.In models that assume a closed-box MDF, the disagree-ments are mostly less than 0.2 dex. Better precisionmetallicity estimates from high N pix pCMDs should bepossible with improvements to the likelihood model (seeSection 5.1).It is important to note that the MIST models (Choiet al. 2016) used in PCMDPy assume a solar α -element ratio, whereas many massive elliptical galaxies (suchas M87 and M49) show strong evidence for super-solar α -abundances and α -abundance gradients (e.g., Sarziet al. 2018). Therefore, the literature metallicities,which include variable abundances, are not perfectlycomparable with our measurements. α -enrichment has been shown to affect broadband col-ors in the RGB and MS in old stellar populations (Dotteret al. 2008), producing redder isochrones. This effect islikely relatively minor, and the offsets could have beenabsorbed in the model fits by the reddening contribu-tion from dust. Because of this, we do not study themeasurements of dust in close detail. DISCUSSION: CAVEATS AND FUTUREAPPLICATIONS5.1.
Likelihood Models for pCMDs
We believe one of the most significant remaining limi-tations to better constrained distances and stellar popu-lations, as well as applying the pCMD technique to moredatasets and with more complex models, is identifyinga better likelihood model for comparing observed andsimulated pCMDs. Comparing distributions of two (or3 log N pix d L i t C1 Single FeH
M87M49NGC3377M31 log N pix d L i t C1 Closed Box MDF
M87M49NGC3377M31
Figure 9 . The measured distances to each galaxy region,compared to the adopted literature values, derived from thesingle metallicity ( top ) and closed box MDF ( bottom ) modelfits and as a function of N pix . The grey shaded region shows a distance uncertainty ( ± . magnitudes). Region C1 inNGC3377 is noted, as its distance measured with the SingleFeH model is inconsistent with the literature value. more) dimensional datapoints is a challenging and un-solved problem, yet this is necessary for comparing dataand model pCMDs.In Cook et al. (2019), we describe our choice to approx-imate the likelihood through comparing the number ofsimulated and data pixels in bins in Hess diagram space,using the size of the bins as a rough control for the ex-pected measurement uncertainty. However, regardlessof whether we assume a Gaussian likelihood model, aPoisson likelihood, or something else (we have experi-mented with many alternative methods), we are left withthe question of whether the likelihood derived in thismanner is statistically meaningful. In other words: canwe reasonably say that one set of parameters are tentimes more likely to have led to the data if this likelihoodmetric is ten times higher than for another set?Even leaving aside stochastic effects (discussed in de-tail in Appendix A of Cook et al. (2019)), we are notconvinced the likelihood framework used here is ideal.Quite often, a very small change in input parameters(such as changing [Fe / H] by 0.02) can result in extremechanges ( (cid:29) × ) to the likelihood. We find this over-specification effect is most significant where the data andmodel are very similar. The likelihood model adoptedhere appears reasonable for discriminating between poorand decent fits, but is oversensitive to small changes toalready reasonable fits. This over specification is primar-ily due to the edge effects introduced by binning, and we suggest that future work be dedicated to identify-ing alternative metrics that do not rely on binning thepCMDs.This behavior led us to the likelihood adjustment pro-cedure described in Section 3.3, which significantly down-weights the samples with the highest likelihood. Philo-sophically, this should be reasonable as we believe thatwe cannot quantitatively trust the relative likelihoods inthe very best-fit points. We believe the samples with thehighest likelihoods are probably somewhat better fits, butwe do not trust the quantitative values. This correctionprocedure will almost certainly result in overestimatingthe uncertainties.A more statistically appropriate likelihood modelwould go a long way towards improving the precisionwith which PCMDPy can constrain distances and stellarpopulations. One possible way forward we have consid-ered is to use a 2-D version of the Kolmogorov-Smirnov(KS) test to evaluate likelihoods. The KS test is a wellestablished metric for comparing fairly arbitrary distri-butions of points in one dimension by comparing thecumulative distribution functions. But the KS test doesnot generalize easily to two or more dimensions, becausethere is no principled way to order points in a cumulativefashion. We have experimented with using the structureof the pCMDs themselves to define an ordering, by cre-ating a dense series of smoothed contours around thedata pCMD points (like those shown in Figures 1 and4), and comparing the relative number of model pointsthat fall within each contour. An alternative likelihoodmodel based on this concept is currently in developmentfor
PCMDPy .5.2.
Applying the pCMD Method at Larger Distances
In this work, we have focused on elliptical galaxieswithin 20 Mpc. As the distance to a galaxy increases,the observed N pix will correspondingly increase, becausethe same number of stars will be confined to fewer pixels.The magnitude of surface-brightness fluctuations scalesroughly as N pix − , due to the Poisson sampling of therare stars. While this implies the precision with whichSFHs and abundances can be inferred will decline atlarger distances, in principle there is no reason why theyshould be completely unmeasurable.To demonstrate this, we also attempted to fit PCMDPy models to archival photometry of NGC 4993, the ellipticalgalaxy at approximately 40 Mpc that was identified asthe host galaxy of GW170817, the first binary neutronstar merger detected in gravitational waves (Cantielloet al. 2018). The data could not be well fit by
PCMDPy , asthe best-fit models predict a distance of less than 10 Mpc,more than × closer than the SBF-measured distance.Despite the failure in our only galaxy studied wellbeyond 20 Mpc, we believe that extending application4 TRGB SBF pCMD (Single) pCMD (MDF)24.224.324.424.524.624.724.824.925.0 D i s t a n c e M o d u l u s ±10% M31
Conn+2016Conn+2012Jensen+2003Tonry+2001 TRGB SBF pCMD (Single) pCMD (MDF)29.629.830.030.230.4 D i s t a n c e M o d u l u s ±10% NGC 3377
Lee+2016bHarris+2007Tonry+2001Ferrarese+2000TRGB SBF pCMD (Single) pCMD (MDF)30.831.031.231.431.631.8 D i s t a n c e M o d u l u s ±10% M87
Lee+2016aBird+2010Blakeslee+2009Mei+2007 TRGB SBF pCMD (Single) pCMD (MDF)30.831.031.231.431.631.832.0 D i s t a n c e M o d u l u s ±10% M49
Mei+2007Blakeslee+2009Jensen+2003
Figure 10 . Comparison of distances to each galaxy by method. Literature distances via the TRGB method and the SBF methodare shown, as well as the weighted average of distances from each region in Figure 9. The pCMD distances show estimatesfrom both assumed metallicity models. For comparison, the gray shaded region shows a ± distance uncertainty aroundthe literature value from Table 1. No TRGB distances to M49 were found in the NED database. All methods are also subjectsystematic calibration uncertainties of order ∼ . magnitudes (not shown). The open circle for NGC 3377 shows the effect ofremoving the apparent outlier, region C1, where we believe the sky background estimate may be contributing to the poor fit. of the pCMD technique to galaxies as far as 50 Mpcshould be possible. The primary practical limitations toextending to larger distances are twofold, and both playa role in the failure of fitting NGC 4993.First, at larger distances, galaxies necessarily extenda smaller area on the sky, reducing the number of adja-cent pixels over which the constant N pix approximationis valid. In NGC 4993, not only is the radial profile aconcern, but we found evidence of notable ( ∼ ) az-imuthal variations in the mean magnitudes of the regionswe extracted, especially in the inner, ∼ kpc, annulus.This is roughly equivalent to averaging over several N pix values, and has the effect of increasing the dispersionin the pCMD. Blanchard et al. (2017) previously de-tected this extended substructure, including radial shellsand strong azimuthal asymmetries, which they argue is evidence of an active history of recent mergers.Secondly, significantly longer exposure times may benecessary to model pCMDs at large distances (and there-fore higher N pix ). Because the average flux per pixelremains constant with distance, so will the error contri-bution from photon noise, such that at fixed exposuretime there will necessarily come a distance beyond whichany intrinsic fluctuations in surface-brightness will bedrowned out by observational errors (see Figure 7 ofCook et al. 2019). In NGC 4993, given the large distanceand relatively short exposure time of the archival pho-tometry, the flux from the galaxy was as little as ofthe estimated sky background in the outer region, andstill only above the background in the inner region.Any uncertainty in the background estimation is there-fore extremely important, and likely also contributed to5 Radius (kpc) [ F e / H ] C1 M87M49NGC3377M31
Figure 11 . Measurements of average [Fe / H] as a function ofradius. Symbols as in Figure 8. There is no evidence for abun-dance gradients, within the large uncertainties. Region C1 ofNGC 3377 is noted, which had a poorly-fit distance. Thehorizontal lines indicate the literature [Z / H] measurementsfrom Table 1. the poor fits.Additional care in the data reduction stage must there-fore be taken when attempting to apply the pCMD tech-nique to more distant galaxies. Highly precise measure-ments of the sky noise will be essential, as will longerexposure times, or stacking of many overlapping expo-sures. Yet, as mentioned in Section 2, stacking manyexposures also requires careful consideration of whichpixels were contributed to by all exposures, in order toaccurately model the properties of photon noise in thesimulations. Nonetheless, these are all issues that canbe addressed with careful attention to the details of thephotometry and proper implementation in the modelingprocess. 5.3. Point-Spread Function Models
In this work, we assume all stars reside at the centersof their respective pixels, while in practice the positionsof stars are randomly distributed within a pixel. Thiseffect can have a noticeable impact on the distributionof pCMDs, as a bright star near the edge of a pixel willdeposit more of its light in neighboring pixels than oneat the center, changing the distribution of the surface-brightness fluctuations.The ideal method to replicate this effect would be tosubsample the simulated pixels (we will call this the subsampled PSF method). In this approach, stars arepopulated (for instance) into a × grid representingsub-pixel positions, their light is distributed using asub-pixel PSF model, and then accumulated into theresulting pixels. This approach adds somewhat to thecomputational cost of simulating a pCMD (effectivelyincreasing N im by a factor of in this example) but would be fairly straightforward.The primary limitation to this approach is access todetailed sub-pixel models for the HST -ACS PSF, as thelatest verison of
Tiny Tim does not provide sub-pixelPSF models appropriate for drizzled
HST photometry.The raw (un-drizzled) photometry is of insufficient qual-ity to perform pCMD analysis due to geometrical dis-tortions, low individual exposure times, and cosmic raycontamination.An alternative approach outlined in Conroy & vanDokkum (2016) and Cook et al. (2019), which we willcall the dithered PSF method, attempts to replicatesub-pixel effects by applying a grid of PSFs, shifted byfractions of a pixel, to different regions of the simulatedimage. This approach therefore assumes all stars in onesection of an image reside at the center of a pixel, whileall stars in another section reside near the edges.To test the applicability of the dithered PSF assump-tion, we compared it to the "correct" subsampled PSFby simulating pCMDs with each technique and assuminga Gaussian PSF, so we can easily generate a sub-pixelPSF model. The dithered PSF approach was found to ar-tificially decrease the variation in simulated pCMDs, bysmoothing the light of stars in some regions of the imageover too many pixels. By contrast, pCMDs generatedusing a single PSF across the entire image (assumingall stars reside at the centers of their pixels) were re-markably similar to those assuming a subsampled PSF,although distinct artifacts in pCMD space appear when N pix (cid:46) . We therefore ignored sub-pixel effects in thiswork, but emphasize that extending the pCMD techniqueto lower N pix (such as in the disk of M31) will require adetailed sub-pixel PSF model for the HST -ACS camera,applicable to drizzled photometry.5.4.
Beyond Intrinsically Smooth Populations
In addition to extending the pCMD technique to largerdistances, there remains the potential of applying it toother nearby systems like spiral galaxies. This workfocused on massive elliptical galaxies because of theirrelative simplicity. The complex structures in disk galax-ies present more significant challenges to the simplifyingassumptions that make fitting
PCMDPy models tractable :namely, that the stars (and dust) in each pixel representindependent draws from the same underlying statisticaldistribution with mean number N pix . These assumptionsare much poorer approximations in late-type galaxies,which have complex features like spiral arms, dust lanes,and star forming clusters.We see two potential approaches that could be usedto apply pCMDs to late-type galaxies. The first wouldinvolve adding more levels of complexity to the forwardmodeling procedure. For instance, each simulated imagecould be represented as having multiple regions with6different N pix (perhaps representing star clusters, spiralarms, or arm gaps) and different dust properties (rep-resenting thick dust lanes). Adding additional modelcomplexity is relatively straightforward: we have alreadyexperimented with allowing a distribution of N pix or atwo-component dust model, and the overhead of creatingmore complex simulated images can be minimized thanksto the GPU-acceleration used in PCMDPy . But addingthese additional components and parameters, many ofwhich are highly degenerate, makes it intractable to fitsuch complex models to data.A better path forward to studying the pCMDs of late-type galaxies may be through more sophisticated datareduction methods , with the goal of extracting relativelyuniform sub-regions, within which the assumptions madeby
PCMDPy are more valid. For instance, the assumptionof constant N pix might be a reasonable approximationof the distribution of stars within a spiral arm gap, orwithin a small region of a spiral arm (if dust lanes andstar clusters are isolated out). Just as we apply PCMDPy to elliptical regions in the galaxies above, disk galaxiescould be partitioned into these sub-regions, perhapsthrough something similar to Voronoi tessellation (e.g.,Cappellari & Copin 2003).5.5.
Complementarity of the PCMD Technique withExisting Methods
In Section 4.2, we demonstrated the capability of thepCMD method to constrain the relative amounts of old( (cid:38) Gyr) star formation in most of the systems studied.The fits were only capable of measuring upper limits onthe younger epochs of star formation, although in a fewcases (primarily in M31) the constraints are sufficient toindicate a significant decrease in the star-formation ratein the most recent ∼ Gyr. Yet we believe it to be afairly generalizeable result that pCMDs, in the currentframework, are limited in their capability of measuringlow levels of young star formation.The primary reason for this is simply the rarity ofyoung stars, on a pixel-by-pixel basis. While massive,young stars are extremely bright and blue, and willtherefore dominate the flux of any pixels they populate,they will contribute to only a very small fraction of thepixels in most pCMDs. Even in the case of a fairly young( τ ∼ Gyr) population, the cumulative contribution ofall stars younger than 10 Myr is around 1 part in .This requires a relatively dense ( N pix > ) region forthere to be one young star per pixel, on average, andeven then the majority will be low-mass, main sequencestars, indistinguishable photometrically from their oldercounterparts. Therefore, only a handful of rare pixelswill be contributed to by young, massive stars, posing achallenge for our current likelihood models to constrainthe amount of young star formation in all but the densest or youngest systems.The pCMD technique, which can constrain old his-tories of star formation, ought to therefore be highlycomplementary to the resolved star and integrated lightmethods, which are useful for constraining young starformation but often insensitive to the oldest ages. Unfor-tunately this complementarity, both in the systems andepochs of star formation that can be studied, makes itchallenging to validate the SFHs measured with thepCMD technique. As first shown in Conroy & vanDokkum (2016), pCMDs can measure old SFHs in thedense bulge of M31, where the resolved star techniquefails to resolve the oldest main sequence turnoff due tocrowding. The entire history of star formation shouldtherefore be recoverable by combining the resolved starand pCMD methods: bright, young stars can be first beresolved and photometered and the remaining crowdedpixels used to study the old star formation with PCMDPy .Integrated light methods are also more sensitive toyoung star formation because they sum the flux overmuch larger regions than individual
HST -ACS pixels,therefore including more of the rare, young stars. Butthese very bright stars typically overwhelm the faintbackground of old main sequence stars (the "outshining"effect, Maraston et al. 2010; Pforr et al. 2012; Sorba &Sawicki 2015), and the ancient SFH is therefore oftendegenerate with other sources of red light, such as dust orAGN activity. The constraints placed on the SFH by theshape of the assumed prior distribution places additionalchallenges on the interpretations (Carnall et al. 2018;Leja et al. 2019). For instance, our measurements of theold SFH in M87 and M49 indicate that the young andold histories of SFHs may not be well represented witha single τ model. If that is true, then integrated lightmethods, which are highly sensitive to the most recentstar formation, may be underestimating the amount ofold SFH when assuming a τ prior. Where high-resolution HST imaging overlaps with multi-band photometry orspectroscopy, combining the pCMD and SED-modelingtechniques could therefore be useful in constraining theentire SFH. SUMMARYIn this work, we present the first application of the pixelcolor-magnitude diagram (pCMD) technique to galaxiesoutside the Local Group. We construct spatially-resolvedpCMDs from archival
HST photometry in three nearbyelliptical galaxies (M87, M49, NGC 3377) and the M31bulge. We align and clean the archival data, and extractthe pCMDS from thin elliptical regions aligned with thegalaxy’s orientation, so that the properties of the datamatch the underlying pCMD modeling assumptions (es-pecially that of fixed N pix ) as closely as possible. We usethe new GPU-accelerated code PCMDPy to fit each dataset7to model pCMDs, and derive Bayesian posteriors overa 5-bin non-parametric star formation history, distancemodulus, and metallicity.We summarize the main results as follows:1. We derive spatially-resolved SFHs in each of thefour galaxies. In the nearest galaxies, M31 andNGC 3377, we constrain the relatively young ( (cid:46) Gyr) history of star formation, and find notableevidence for radial gradients in the SFH. The starsin the inner regions formed significantly earlier thanin the outer regions. In M49 and M87, only thevery oldest ages of star formation are constrained.2. At similar radii, the bulge of M31 has a youngerSFH than NGC 3377, which in turn shows evidenceof more recent star formation than M87 or M49.3. We recover distance estimates to the four galaxies,in good ( ∼ ) agreement with the tip of thered-giant branch (TRGB) and surface-brightnessfluctuation (SBF) techniques through combiningthe individual distance estimates in each region.4. pCMDs are more sensitive to the oldest epochs of star formation, and so the pCMD techniqueshould be highly complementary to existing SFHmeasurement techniques that primarily constrainyoung star formation.We thank Peter Blanchard and Ashley Villar for help-ful discussions about the details of HST photometryand data reduction. B.C. acknowledges support fromthe NSF Graduate Research Fellowship Program undergrant DGE-1144152. This work is supported in part byHST-AR-14557. The computations in this paper wererun on the Odyssey cluster supported by the FAS Divi-sion of Science, Research Computing Group at HarvardUniversity.
Software:
This research has made use of NASA’sAstrophysics Data System, as well as the following soft-ware packages: PyCUDA (Klöckner et al. 2012), DynestySpeagle (2019), NumPy (van der Walt et al. 2011), Mat-plotlib (Hunter 2007), IPython (Pérez & Granger 2007),Jupyter (Kluyver et al. 2016), SciPy (Jones et al. 2001),Pandas (McKinney 2010), and Astropy (The Astropy Col-laboration et al. 2013, 2018).REFERENCES
Barbary, K., & Kyle. 2016, The Journal of Open Source Software,1, 58Bertin, E., & Arnouts, S. 1996, Astronomy and AstrophysicsSupplement Series, 117, 393Binney, J., & Merrifield, M. 1998, Galactic astronomy (PrincetonUniversity Press), 796Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, TheAstrophysical Journal, 694, 556Blanchard, P. K., Berger, E., Fong, W., et al. 2017, TheAstrophysical Journal, 848, L22Bothun, G. D. 1986, The Astronomical Journal, 91, 507Cantiello, M., Jensen, J. B., Blakeslee, J. P., et al. 2018, TheAstrophysical Journal, 854, L31Cappellari, M., & Copin, Y. 2003, Monthly Notice of the RoyalAstronomical Society, Volume 342, Issue 2, pp. 345-354., 342,345Carnall, A. C., Leja, J., Johnson, B. D., et al. 2018, TheAstrophysical Journal, 873, 44Choi, J., Dotter, A., Conroy, C., et al. 2016, The AstrophysicalJournal, 823, 102Conn, A. R., McMonigal, B., Bate, N. F., et al. 2016, MonthlyNotices of the Royal Astronomical Society, 458, 3282Conroy, C., Dutton, A. A., Graves, G. J., Mendel, J. T., & vanDokkum, P. G. 2013, The Astrophysical Journal, 776, L26Conroy, C., & van Dokkum, P. 2012, The Astrophysical Journal,760, 71Conroy, C., & van Dokkum, P. G. 2016, The AstrophysicalJournal, 827, 9Cook, B. A., Conroy, C., Dokkum, P. v., & Speagle, J. S. 2019,The Astrophysical Journal, 876, 78Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, TheAstrophysical Journal Supplement Series, 200, 18Dolphin, A. E. 2002, Monthly Notices of the Royal AstronomicalSociety, 332, 91 Dotter, A., Chaboyer, B., Jevremović, D., et al. 2008, TheAstrophysical Journal Supplement Series, 178, 89Ferrarese, L., Mould, J. R., Kennicutt, Jr., R. C., et al. 2000, TheAstrophysical Journal, 529, 745Fruchter, A., Hack, W., Dencheva, N., Droettboom, M., &Greenfield, P. 2010, in STSCI Calibration WorkshopProceedings, ed. S. E. Deustua & C. Oliveira, Baltimore, MD,376Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, TheDrizzlePac Handbook (Baltimore, STScI: Space TelescopeScience Institute)Holmberg, J., Nordström, B., & Andersen, J. 2007, Astronomy &Astrophysics, 475, 519Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MonthlyNotices of the Royal Astronomical Society, 445, 581Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90Jee, M. J., Blakeslee, J. P., Sirianni, M., et al. 2007, Publicationsof the Astronomical Society of the Pacific, 119, 1403Jones, E., Oliphant, T., Peterson, P., & Others. 2001, SciPy:Open source scientific tools for PythonJoye, W., & Mandel, E. 2003, in Astronomical Data AnalysisSoftware and Systems XII, Vol. 295 (Astronomical Society ofthe Pacific (ASP)), 489Klöckner, A., Pinto, N., Lee, Y., et al. 2012, Parallel Computing,38, 157Kluyver, T., Ragan-Kelley, B., P{\’e}rez, F., et al. 2016, inPositioning and Power in Academic Publishing: Players, Agentsand Agendas, 87Krist, J. E., Hook, R. N., & Stoehr, F. 2011, in Proceedings of theSPIE, Volume 8127, id. 81270J (2011)., ed. M. A. Kahan, Vol.8127, 81270JLanyon-Foster, M. M., Conselice, C. J., Merrifield, M. R., et al.2007, MNRAS, 380, 5718