The Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER) II. The Spatially Resolved Recent Star Formation History of M33

We measure the spatially resolved recent star formation history (SFH) of M33 using optical images taken with the Hubble Space Telescope as part of the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER) survey. The area covered by the observations used in this analysis covers a de-projected area of $\sim$38 kpc$^{2}$ and extends to $\sim$3.5 and $\sim$2 kpc from the center of M33 along the major and semi-major axes, respectively. We divide the PHATTER optical survey into 2005 regions that measure 24 arcsec, $\sim$100 pc, on a side and fit color magnitude diagrams for each region individually to measure the spatially resolved SFH of M33 within the PHATTER footprint. There are significant fluctuations in the SFH on small spatial scales and also galaxy-wide scales that we measure back to about 630 Myr ago. We observe a more flocculent spiral structure in stellar populations younger than about 80 Myr, while the structure of the older stellar populations is dominated by two spiral arms. We also observe a bar in the center of M33, which dominates at ages older than about 80 Myr. Finally, we find that the mean star formation rate (SFR) over the last 100 Myr within the PHATTER footprint is 0.32$\pm$0.02 M$_{\odot}$ yr$^{-1}$. We measure a current SFR (over the last 10 Myr) of 0.20$\pm$0.03 M$_{\odot}$ yr$^{-1}$. This SFR is slightly higher than previous measurements from broadband estimates, when scaled to account for the fraction of the D25 area covered by the PHATTER survey footprint.


INTRODUCTION
Corresponding author: Margaret Lazzarini mlazz@caltech.edu The recent star formation history (SFH) of a galaxy allows us to look back in time to see the effects of the astrophysical phenomena that have influenced its evolution. If this recent SFH can also be spatially resolved, we can see the impact of the astrophysics that controls the stellar structure of the galaxy including the forma-tion and evolution of spiral arms and the connection between local star formation events and the surrounding dust and gas from which new stars form.
There are several methods for deriving the recent SFH of a galaxy. For studies of distant galaxies, integrated light from unresolved stars can relate a galaxy's current star formation rate to other galaxy-scale properties such as its mass, luminosity, and color. While these studies are useful for obtaining statistically significant results with large samples of galaxies, they do not allow us to study the physics surrounding star formation within each galaxy in detail. For this, we must turn to the more local universe where we can resolve star formation across individual galaxies rather than approximating their physical properties with a handful of parameters.
Most methods for measuring the recent SFH of local galaxies rely on star formation rate (SFR) indicators, which trace the presence of young, newly formed stars through a variety of signatures. These tracers can be both direct and indirect, each with a characteristic time scale. Ultraviolet (UV) emission from young massive stars traces star formation within the last 100−200 Myr, corresponding to the lifetimes of these stars (Kennicutt & Evans 2012). As a direct tracer of star formation, UV emission has been widely used to measure star formation rates in nearby galaxies, especially after the launch of GALEX (e.g., Gil de Paz et al. 2009). However, UV flux is attenuated by interstellar dust, making it more of a lower limit on the SFR (e.g., Bell 2003). Hα emission from gas ionized by young O-stars traces very recent star formation within the last ∼5 Myr (Byler et al. 2017). Hα is thus an indirect tracer of young, massive stars (>15 M ), but is one of the best available measure of the instantaneous SFR in external galaxies. Dust attenuation (Kennicutt & Evans 2012) and the leakage of ionizing photons (Lee et al. 2009;Choi et al. 2020) are sources of systematic uncertainty in Hα measurements.
A more direct method for measuring the star formation rate is to not just confirm the presence of young, massive stars, but instead to actually count the number of stars of various masses directly and use this distribution to infer time-resolved star formation rates that would produce the observed population (e.g., Dolphin 2002;Weisz et al. 2008;Williams et al. 2011;Kennicutt & Evans 2012). This method requires deep, high resolution observations to resolve individual stars for a full population, but provides the information necessary to produce a spatially-resolved star formation history. We can use resolved stars to construct color-magnitude diagrams (CMDs) for specific spatial regions within galaxies. The distribution of stars on the CMD is a prod-uct of the galaxy's history of star formation, the evolution of the galaxy's metallicity, and the initial mass function (IMF), all observed through the distribution of dust. Like all SFH measurement techniques, modelling the CMD relies on assumptions about the IMF, stellar evolution and atmospheric models, binary fraction and dust. While CMD fitting is is not the only method we currently have to produce time-resolved SFR measurements -SED-fitting can also produce time-resolved SFR measurements -it does produce the highest timeresolution in these measurements.
Much of the early work using CMD fitting on galaxies in the local universe focused on low mass galaxies and smaller portions of large galaxies. Dwarf galaxies are an excellent target for resolved stellar populations work because of their abundance in the Local Volume, making them close enough to resolve their individual stars (e.g., Harris & Zaritsky 2004, 2009McQuinn et al. 2010;Weisz et al. 2011Weisz et al. , 2014Geha et al. 2015;Cignoni et al. 2019). Studies focused on larger galaxies have often relied on small fields across the disk to answer questions about galaxy evolution, such as inside-out growth. However, these small fields make resolving the temporal evolution of large scale structures such as spiral arms or bars difficult (e.g., Williams 2002Williams , 2003Brown et al. 2006Brown et al. , 2007Brown et al. , 2008Williams et al. 2009Bernard et al. 2012Bernard et al. , 2015Choi et al. 2015).
CMD modelling has successfully derived spatially resolved SFHs for well-resolved galaxies with wider field coverage with multiple HST fields. There have been several significant applications of this SFH measurement technique in nearby galaxies including M31 (Lewis et al. 2015;Williams et al. 2017), M33 (Williams et al. 2009), M81 (Choi et al. 2015), NGC 300 (Gogarten et al. 2010), and NGC 2403 (Williams et al. 2013). The common distance of stars in these nearby galaxies allows for a more simplified CMD modelling approach over stellar populations in the Milky Way. Lewis et al. (2015) provided the first spatially resolved star formation history of a significant portion of an L * galaxy, M31, using photometry from the Panchromatic Hubble Andromeda Treasury (PHAT). Their fits were optimized for main sequence stars, providing a recent star formation history for roughly 9000 regions measuring ∼100 pc by 100 pc across a 0.5 square degree area in the northern third of the star forming disk within the last ∼500 Myr. Williams et al. (2017) provided an analogous measurement of the spatially resolved ancient star formation history of M31 by focusing their fits on older stellar populations. Both studies capitalized on the rich dataset from PHAT, which provided six-band photometry from the near-infrared through near-UV for over 100 million individual stars in M31 Williams et al. 2014).
The spatially resolved SFH maps of M31 derived via CMD-modelling have had several significant applications in subsequent work. The Lewis et al. (2015) maps were important in establishing that M31 had a disk wide burst of star formation when it merged (D'Souza & Bell 2018;Hammer et al. 2018). By comparing the spatially resolved SFH maps derived via CMD-modelling with SFRs derived from various indicators, Lewis et al. (2017) found that the SFR indicators in M31 were under-estimating the total amount of star formation. Díaz-Rodríguez et al. (2018) used SFH maps to measure the progenitor mass distribution for core collapse supernovae in M31 and M33 using the SFH maps of Lewis et al. (2015) and Jennings et al. (2014). The Lewis et al. (2015) SFH maps in M31 have also been used to measure the age distribution of high mass X-ray binaries in M31 Lazzarini et al. 2018Lazzarini et al. , 2021.
The insight into M31's evolution gained from the PHAT survey has been incredibly valuable. In this paper, we derive the spatially resolved recent SFH of M33. Though the techniques are the same, M33 and M31 are very different galaxies (e.g., in mass, structure, and starforming activity), providing important context to the results from PHAT. The star formation rate intensity of M33 is higher than M31 (Verley et al. 2009;Lewis et al. 2015), allowing us to study high intensity star formation regions. M33 also has well known stellar and gas-phase metallicity gradients (Cioni 2009;Magrini et al. 2009Magrini et al. , 2010Beasley et al. 2015;Toribio San Cipriano et al. 2016;Lin et al. 2017), which are similar to the LMC (Carrera et al. 2008). Unlike M31, M33 has likely not undergone a recent significant merger event, though recent spectroscopic observations may hint at a centrallyconcentrated population of accreted stars (Gilbert et al. 2021). Its warped outer disk does suggest a mild encounter with M31 in the last 1-3 Gyr (Putman et al. 2009), but proper motions measurements indicate that M33 is likely on its fits infall into the M31 system (e.g., Patel et al. 2017b,a;van der Marel et al. 2019). This makes it an excellent target for CMD-fitting based SFH measurements, allowing results to be compared directly to the SFHs of disrupting Magellanic spirals such as the LMC (Mazzi et al. 2021), and undisturbed field spirals such as NGC 300 (Gogarten et al. 2010).
M33 is of comparable mass and metallicity to the Large Magellanic Cloud (LMC; M ∼ 3×10 9 M , [M/H] ∼ −0.5; e.g., van der Marel et al. 2002;Carrera et al. 2008;van der Marel et al. 2012;Beasley et al. 2015), making it one of the most massive satellite galaxies in the Local Group. It does not have a significant bulge component (Kormendy & McClure 1993;McLean & Liu 1996) and has been observed to host a weak bar (Elmegreen et al. 1992;Regan & Vogel 1994;Block et al. 2004;Hernández-López et al. 2009;Corbelli & Walterbos 2007), both of which reduce confusion between bulge and disk stellar populations. Its proximity allows us to resolve stars down to the ancient main sequence (Williams et al. 2009). In fact, we can resolve its stellar populations better than we can in M31 due to its similar distance but lower stellar surface density.
There has been significant previous work to measure the spatially-resolved SFR of M33 using star formation indicators across the electromagnetic spectrum and resolved stellar populations, and we briefly review them here. These include the spatially resolved SFR of M33 using broadband flux SF indicators and resolved stellar populations CMD-modelling.
The spatially resolved SFR of M33 has been previously studied using star formation indicators including FUV, Hα, and 24 µm emission. Verley et al. (2009) measured the global SFR of M33 over the last 100 Myr using extinction corrected FUV and Hα flux and detected a radial decline in SFR. Verley et al. (2007) measured the global and local SFR of M33 using 24 µm emission from HII regions. Far-infrared emission from dust is an important SF indicator because much of the FUV flux from young stars is absorbed and re-emitted by surrounding dust (Boquien et al. 2010(Boquien et al. , 2011. There can be large discrepancies between SF measurements with FUV, Hα, and far-infrared SF indicators (Boquien et al. 2015). This highlights the importance of CMD modelling, which does not rely on the calibration of SF indicators and can also produce time-resolved measurements. However, calibrating SF indicators with nearby galaxies such as M33 is very important work because we can only measure SFHs with CMD modelling for galaxies that are near enough to resolve their stellar populations.
CMD-fitting has previously been used to measure the recent SFH of M33 in targeted regions of the disk, confirming a mean age and metallicity gradient in the disk, suggesting inside-out growth (Barker et al. 2007;Williams et al. 2009;Barker et al. 2011). CMD-based SFH measurements of SFHs around supernova remnants in M33 and M31 suggest that some massive stars may not explode as supernovae (Jennings et al. 2014). These CMD-based SFH measurements were targeted and did not cover the full disk, or a large contiguous area of the disk of M33. The PHATTER data now allow such a measurement for the first time, opening the door to studies of temporal evolution in morphological features.
In this work we generate maps of SFH within the last ∼630 Myr in M33 using optical photometry from the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER) survey . Our approach mirrors that in the Lewis et al. (2015) analysis with the PHAT data in M31, allowing the two measurements to be directly compared. This work is one of the first projects using the rich PHATTER dataset to study the star formation, dust, morphology, and more with photometry for over 22 million individual stars in M33.
In this paper we use optical photometry from the PHATTER survey  to derive the spatially resolved recent star formation history of M33. In Section 2 we present the optical photometry used as well as our method for creating artificial stars to measure photometric errors. In Section 3 we present our method for deriving the spatially resolved recent SFH of M33. In Section 4 we present the resulting SFH maps and in Section 5 we discuss our results. In Section 6 we summarize our findings.
For the analysis presented in this paper, we assume a distance modulus for M33 of 24.67 (de Grijs & Bono 2014), which is equivalent to a distance of ∼859 kpc. We use the measured D25 ellipse presented by Gil de Paz et al. (2009) which is centered at 1 h 33 m 50 s .9, 30 • 39 35 .8. We adopt a position angle (P.A.) of 201 • and an inclination angle of 55 • (Koch et al. 2018).

PHATTER DATA
We derive the spatially-resolved recent star formation history of M33 using optical photometry from the PHATTER survey 1 . The optical-only PHATTER survey covers a projected area of ∼ 300 arcmin 2 , or ∼38 kpc 2 de-projected, in the center of the star forming disk of M33, as shown in Figure 1. While the six-band PHATTER photometry catalog contains over 22 million individual stars and includes six filter photometry from the near-infrared to near-ultraviolet, as discussed in detail in Williams et al. (2021) the data we use in our analysis comes from the ACS (optical; F475W and F814W) data alone. This catalog from the optical-only imaging data contains ∼25 million stars total with ∼15 million of these having reliable measurements (i.e., stars that meet our stringent "gst" quality cuts which are described in Section 2.1) in both the F475W and F814W bands. The use of the optical-only catalog was critical for creating the artificial star tests (ASTs) in a reasonable time and creating a catalog of homogeneous exposure times across the survey footprint. We compute the SFH for all of the 1 PHATTER Data, DOI: 10.17909/t9-ksyp-na40 area within the PHATTER survey outline indicated in Figure 1.

Photometry and Catalogs
We use optical photometry in the F475W and F814W bands for our measurements. We use the "gst" catalogs, which contain "good stars" that pass cuts related to their signal to noise, crowding, and sharpness as described in Williams et al. (2021). Specifically, all stars in the "gst" catalog have S/N>4, sharpness 2 <0.2, and crowding<2.25 in both bands. These cuts are designed to leave only the stars with the highest quality photometry for our measurements. For our analysis, we divided the "gst" photometric catalog into 24 square regions (roughly 100 pc on a side), along lines of constant R.A. and Dec., producing 2005 regions across the survey area. The mean number of "gst" stars in our analysis regions is ∼4000, and the mean number of "st" stars (stars with S/N≥4) is ∼13,000.

Artificial Star Tests
It can be challenging to accurately measure the magnitudes of stars in crowded regions. Faint stars can blend together, biasing the measured flux from other stars within the same region, such that magnitudes of stars in regions of high stellar density will have more biased, noisier measured flux. We quantify the degree of bias, noise, and completeness in our photometric measurements using ASTs.
We repeatedly insert artificial stars into an image and then run the image through our standard photometric pipeline, as described in Williams et al. (2021). The measured difference between the input and recovered magnitude of these stars quantifies the effect of crowding in that region. The same test can be used to set the completeness. We calculate the magnitude at which 50% of the artificial stars are recovered and use this as our limiting magnitude when we measure the SFH in each region. Because these effects are sensitive to both stellar flux and stellar density, we want to repeat this process for stars across the CMD for a range of local stellar densities.
Rather than fully populating each pixel with the required number of artificial stars, we take advantage of the fact that the photometric effects of crowding are primarily density-dependent. We can combine ASTs from across the galaxy, as long as the background stellar density is similar. We quantify local density using the number of stars in the "st" catalog (stars with S/N≥4) per square arcsecond. We use the "st" rather than "gst" catalog to better account for the real stellar density in the region, not just the density of stars that meet our stringent "gst" criteria. We grouped artificial star tests by The PHATTER survey footprint overlaid on a color optical image of M33 (generated with R, V, and B band images from the Local Group Galaxy Survey; Massey et al. 2006). The outline was generated with the ACS exposure map from the PHATTER survey and includes all regions for which star formation histories were calculated for this paper. The full survey outline was divided into 2005 individual regions that measure 24 on a side, using a grid of constant R.A. and Dec. right: Image showing the ACS exposure map from the PHATTER survey. The exposure map matches the area of the PHATTER survey outline shown in the left panel of this figure. The exposure map is color coded by total exposure time, measured in ks. stellar density into 16 stellar density bins to increase the number of tests that can be used for each of the SFH regions. For each of the ∼2000 regions in our SFH analysis, we used a set of ASTs that were generated for regions of similar stellar density. The number of ASTs and the 50% completeness limits for each stellar density bin in Table 1. The 50% completeness limits correspond to an old main sequence turn off (MSTO) age of about log(age/yr)=9.5 at the faint end and about log(age/yr)=9.2 at the bright end. This suggests that our recent ( 630 Myr) SFH maps should be complete, even in the most crowded regions.

DERIVATION OF THE SFHS
We use the CMD-fitting code MATCH to derive the SFHs (Dolphin 2002). MATCH finds the combination of single-burst stellar populations that best reproduces the observed CMD. MATCH incorporates photometric errors and completeness that are measured in the ASTs (see Section 2.2) and fits for dust extinction. The best fit SFH is then determined using a Poisson maximum likelihood technique. For more details, see Dolphin (2002) and Lewis et al. (2015).
Because our focus is on the recent SFH only, we exclude regions of the CMD occupied by older, redder stars. The atmospheres and evolution of main sequence stars are generally better understood than later evolutionary stages, making the SFH recovery for young ages reasonably secure. We exclude the area of the CMD dominated by these older ages, such as the vast majority of red giant branch and red clump stars, with both F475W magnitude greater than 21 and F475W-F814W color greater than 1.25.
Our choice of exclusion region purposefully includes main sequence stars and young red and blue He burners in our fits. Red and blue He burners are useful for age dating at young ages (e.g., McQuinn et al. 2012). We chose the blue edge of our exclusion region (at F475W−F814W=1.25) with the goal of excluding red giant branch and red clump stars, without also excluding the most reddened young stars. If we moved our exclusion region too far towards the blue we would risk excluding the young, reddened population from our fits, which would potentially bias the results towards recovering star formation in the affected age ranges only in regions with lower-extinction values. The average number of stars we fit on the CMD for a given analysis region is ∼4000, with a standard deviation of ∼1500 stars. The MATCH fitting of a given region is essentially finding the intrinsic distribution of young stars and dust that reproduces the luminosity function and color distribution of the main sequence, given photometric selection, noise, and biases from the ASTs.
In our fits we define the input the age bins, metallicity, distance, extinction, binary fraction, and choice of initial mass function (IMF). Where possible, we adopt the same model parameters that Lewis et al. (2015) used to fit the SFHs in M31 to make our results directly comparable. We set a fixed distance modulus of 24.67 (de Grijs & Bono 2014), use a binary fraction of 0.35, and draw the mass of the secondary from a uniform distribution and Kroupa IMF (Kroupa 2001). We use the same set of age bins as Lewis et al. (2015), which span from 6.6 to 10.15 log(years) with a spacing of 0.1 dex, with the last bin spanning from log(years)=9.9−10.15. The youngest age of the isochrones used in this work is 4 Myr. We calculate the stellar mass formed in the youngest time bin (between 6.6 and 6.7 log(years) ) for each spatial region, and then derive the SFR that extends to the present day by dividing that stellar mass by the total time between the upper bin edge (log(years)=6.7) and the present day (0 Myr). The fine time resolution at young ages, within the last 25 Myr or so, results in high uncertainties in our SFR measurements, due to the small number of stars of these ages on the CMDs and theoretical uncertainties. Users of the SFH measurements presented in this paper may combine time bins in cases where they would prefer lower resolution with smaller errors.
We adopt metallicity constraints for the stellar model set used. Although main sequence evolutionary models are better constrained than other phases, there are differences between model sets, mostly at high mass. We perform our calculations using both the Padova (Marigo et al. 2008)  To generate our SFHs, we must generate model CMDs with stars of various ages and metallicities distributed over a two dimensional parameter space of age and metallicity. While age-metallicity degeneracy is not as strong for main sequence stars, the metallicity will strongly effect the color of the main sequence (e.g., Gallart et al. 2005). M33 has a well-established stellar and gas-phase metallicity gradient (Cioni 2009;Magrini et al. 2009Magrini et al. , 2010Toribio San Cipriano et al. 2016;Lin et al. 2017) and we explore the effects of metallicity on our results in Section 4.4.

Choice of Region Size
To generate the spatially resolved recent star formation history of M33, we use a region size of ∼ 24 that is approximately 100 pc on a side. This region size was chosen to match the size of regions that were used to measure the spatially resolved recent star formation history of M31 to make our results directly comparable. As described by Lewis et al. (2015), the 100 pc region size was chosen because it bridges the gap between previous work on Galactic pc scale star formation and star formation in more distant galaxies on kpc scales. In addition to scientific interest, the region size was chosen to provide a similar number of stars for the CMD for each individual region and to fit within the computational constraints of our analysis. Given the comparable distance of M33 to M31 and the analogous survey design, we chose to use the same region size for our analysis, which also eases comparison between the M33 and M31 results.
We generate our grid of analysis regions along lines of constant R.A. and Dec. Because of this, along the edges of the PHATTER survey footprint, some of our regions have partial coverage in the optical-only catalogs. We calculate the area of each SFH region covered by the optical-only PHATTER survey and include that value in units of [arcsec 2 ] in Tables 2 and 3. We use the de-projected area of each SFH region covered by the optical-only PHATTER catalog to generate figures where the SFR is shown in units of M yr −1 kpc −2 .

Extinction
M33 is lower metallicity and more face on (lower inclination) than M31 and thus has only moderate dust extinction (Corbelli 2003). It does, however, have dust that needs to be accounted for when modeling the CMD. We include dust in our CMD fitting by including two parameters that describe the dust in each SFH region: foreground extinction (A V ), which reddens all the stars uniformly, and differential extinction in star forming regions (dA V ). We employ the same method used by Lewis et al. (2015) to recover extinction in the spatially resolved recent SFH measurement using the PHAT data in M31 in an effort to make our measurements directly comparable. This extinction model is adequate for young stars (age<1 Gyr), which are confined to the plane and tend to experience a roughly tophat differential reddening distribution (e.g., Dolphin et al. 2003;Weisz et al. 2014). Older stellar populations require more complicated fitting for extinction which follows a log-normal distribution Williams et al. 2017). In MATCH, extinction is applied to all stars in a uniform distribution between A V and A V + dA V . Thus, A V shifts the entire CMD fainter and red-ward, whereas dA V broadens and dims the main sequence. We briefly describe the method we use for fitting A V and dA V below.
We determine the best values for A V and dA V by running our SFH calculations multiple times over a grid of A V and dA V . We limit our A V values to a range of [0.0,1.5] and our dA V values to a range of [0.0,2.5]. First, we run the SFH calculations over a grid with an initial spacing of 0.2 in both values. We calculate the likelihood of each (A V ,dA V ) pair and then choose the pair with the highest likelihood value. Then we create a finer grid with 0.05 spacing in both values within 2σ of the best (A V ,dA V ) pair and repeat our maximum likelihood technique to identify the best fit (A V ,dA V ) for the region. This method is described in more detail in Appendix B of Lewis et al. (2015).

Uncertainties
Uncertainties in our measurement of the SFH of M33 come from a combination of random and systematic uncertainties, and dust.
The random uncertainties were calculated using a hybrid Monte Carlo (MC) process (Duane et al. 1987) implemented as part of the MATCH software package (Dolphin 2013). These random uncertainties scale with the number of stars on the CMD for a given region, with more sparsely populated CMDs resulting in higher random uncertainties. Rather than generating SFHs on a grid like we do for the A V , dA V calculations, the hybridMC routine uses a Monte Carlo routine that explores high dimensional parameter spaces more efficiently than a traditional Metropolis-Hastings Monte Carlo routine. The hybrid MC routine generates 10,000 possible SFHs for the CMD from a given analysis region and the density of the generated SFHs in parameter-space is proportional to the probability density. The errors on the SFHs can then be measured by identifying the boundaries of the region containing 68% of the samples. This method accounts for cases where the maximum likelihood SFH has zero star formation in individual time bins, allowing the uncertainties of the SFR on those time bins to be estimated.
To assess systematic uncertainties due to our choice of stellar models, we ran our full set of SFH calculations with two sets of stellar models: the Padova models with complete AGB tracks (Marigo et al. 2008;Girardi et al. 2010) and the MIST stellar models (Dotter 2016). Differences between stellar models mainly stem from how they model different aspects of stellar evolution including mass loss, convection, and rotation (Conroy 2013). These differences tend to be greater in post-main sequence phases of stellar evolution and for the most massive stars.
We consider uncertainties due to dust content in each individual SFH region. We expect to underestimate the star formation rate at young ages in some spatial regions due to embedded star formation. We do not expect this underestimation to affect all spatial regions equally, due to the uneven distribution of dust across the PHATTER survey footprint. We address this uncertainty in more detail and calculate the fraction of star formation we expect to have missed due to embedded star formation in Section 4.3.
Lastly, we acknowledge that our choice of binary fraction could be an additional source of systematic uncertainty in our SFH calculations. For consistency, we adopt the same binary fraction (0.35) used by Lewis et al. (2015) in deriving the recent SFH for M31 using the PHAT data. As shown in Lewis et al. (2015), uncertainties introduced by choice in binary model are small compared to dust variation uncertainties.

Reliable Age Range of SFH Measurements
Our SFH measurements are not sensitive to the older stellar populations that were excluded from the red side of the CMD in our fits. Our SFHs are therefore at best only reliable back to the ages of the oldest main sequence stars analyzed in each region. However, the true age back to which our fits are reliable can vary by region due to the effects of crowding and dust extinction.
To determine the sensitivity of our SFH fits as a function of lookback time, we created artificial CMDs with MATCH, which we present in Figure 2. We created artificial CMDs with constant SFR across age bins and solar metallicity and varied the stellar density and dust extinction, which are the two main factors that would affect the reliability of our SFHs with respect to lookback time. We created artificial CMDs across a grid of stellar density and dust extinction. In the first three rows of Figure 2 we show artificial CMDs for low, medium, and high stellar density regions. We repeat the artificial CMDs for the high stellar density region in the fourth row. The points on the artificial CMDs are color-coded by stellar age. All stars older than 800 Myr are plotted in black and stars younger than 800 Myr are color coded by age according to the colorbar on the right-hand side of the figure.
The low, medium, and high stellar density regions were created with stellar densities of 0−10, 22−24, and 34−38 stars per square arcsecond, respectively. To generate artificial CMDs at each stellar density, we used the completeness limits measured for each stellar density bin and used ASTs for the given stellar density. Our artificial CMDs were created with three levels of dust extinction: no dust extinction ( . Median and high dust extinction are defined from the range of A V and dA V values we measured in our SFHs.
In the first three rows of Figure 2, we plot 630 Myr isochrones reddened with the A V specified for each column with a thick black line. We also include a thinner, dotted line showing the same age isochrone reddened with A V + dA V , to show the most extreme reddening applied to stars in our fits. In the fourth row, we plot a 630 Myr isochrone reddened with the A V specified for each column and a 1 Gyr isochrone with the same reddening applied.
We find that at low stellar density and extinction our photometry depth allows us to observe main sequence stars back to about 1 Gyr. However, when the depth is limited by crowding in the high stellar density regions, the main sequence turn-off for stars older than ∼800 Myr falls at or near the magnitude limit, meaning that our SFH measurements are not reliable at these ages. We thus adopt a conservative age limit of 630 Myr for our scientific analysis in this paper, although lower density and less-extincted regions may have reliability back to older ages. This reliability limit is older than the one adopted by Lewis et al. (2015) in the spatially resolved recent SFH of M31 (400 Myr). This difference is caused by the higher crowding of M31, due to its more edge-on inclination. Higher crowding results in brighter magnitude limits, resulting in SFH measurements that are reliable back to younger ages.
We include figures showing CMDs with our data for three regions (low, medium, and high stellar density) and compare these CMDs with model CMDs in Appendix A. , shown in the left, middle, and right column, respectively. Each row shows CMDs with different stellar densities, specified in the upper left corner of the plots in the first column. For more details in the exact densities used, see Section 3.4. The points on each CMD are color-coded by stellar age in Myr. All stars older than 800 Myr are plotted in black. The gray box outlines the range of colors and magnitudes excluded in our fits. In the first three rows we plot two 630 Myr isochrones reddened with the given column's AV (thick solid line) and AV + dAV (thin dotted line). In the fourth row, we plot artificial CMDs for high stellar density and 630 Myr and 1 Gyr isochrones, both reddened with AV . For further discussion, see Section 3.4.
We ran tests where known SFH and extinction were used to create synthetic CMDs, from which the SFH and extinction were fit. We ran these tests with realistic SFHs and with single bursts of star formation. We discuss these tests and their implications in detail in Appendix B.

RESULTS
In this section we present the integrated and spatially resolved SFHs for M33. We present the full SFH results for each spatial region in Tables 2 and 3, where we present column descriptions and the full SFR measurements and the coordinates of the axes of each individual region. We include the SFH for the first six spatial regions in Table 3 and include the full table for all spatial regions in the digital version of this paper. Unless otherwise specified, figures and tables showing SFH results were created with the Padova model measurements.
In Figure 3, we explore the effect of stellar density on our SFH measurements for individual regions. In panel (a) we include a color image of M33 with the PHATTER survey footprint outlined and the location of three of our analysis regions outlined and labeled with their region id numbers. These regions were chosen for their varying stellar densities. Region 1775 is a lower stellar density region (∼7 stars arcsec −2 ), region 1165 is a medium stellar density region (∼24 stars arcsec −2 ), and region 1095 is a high stellar density region (∼34 stars arcsec −2 ). For each region, we include a CMD of the "gst" stars in that region, with the region of the CMD excluded from our fits plotted as a gray box. Stars that lie outside of the excluded region were used to measure the SFH in that region. In panels (b, c, d) we include a color composite image created with the PHATTER survey mosaics (red=F160W, green=F814W, blue=F475W). The SFH plotted for each region includes a histogram showing the best-fit SFR for each time bin, with random uncertainties. The uncertainties presented for each region are typical of those for regions of similar stellar density. While the increased stellar density can affect our completeness limits, the increased number of stars being fit on the CMD results in smaller error bars on our SFR measurements.
The co-added SFH measurements for all sub-regions reveal the integrated recent SFH of M33 within the PHATTER footprint. We add the stellar mass formed within all spatial bins for the same time bin. The integrated SFH within the PHATTER footprint obtained with the Padova models is listed in Table 4. We list the start and end of each time bin in both log(years) and Myr and include uncertainties, which include both random and dust uncertainties. We plot the spatially in-tegrated SFH in Figure 4, where we present the best-fit SFH obtained with both the Padova and MIST models as a function of time. The dark histogram represents the Padova models and the lighter pink histogram represents the MIST models.
The horizontal dashed lines in the lower panel of Figure 4 (which overlap) indicate the mean SFR within the last 630 Myr (∼ 0.34 M yr −1 ) derived for the PHATTER survey area for the MIST and Padova models. While their SFR measurements within the youngest time bins differ, the mean SFR within the last 630 Myr measured with the two models differs by <1%.
The MIST and Padova models agree within their random uncertainties for most time bins older than about 20 Myr, indicating that random uncertainties dominate the errors on the SFR at these ages. The largest systematic differences in our SFR measurements are in the youngest time bins, where models for the most massive stars generate most of the uncertainty. At the youngest ages (<15 Myr), the model sets differ by up to a factor of two in the measured SFR.
We also present the SFH (for the Padova models) plotted linearly in Figure 5, to better represent the resolution of the youngest time bins. We find that the SFR of M33 is fairly constant over the past ∼500 Myr. Both the Padova and MIST model measurements suggest some recent elevation in the SFR, but they place it at slightly different times and amplitudes. The Padova models find two peaks in the SFR within the past 100 Myr, around 5−10 Myr and 70−100 Myr ago, with an elevation of a few times 10 −2 M yr −1 above the mean SFR over the past 630 Myr.

SFH Maps
In Figure 6, we show the spatially-resolved SFH, plotting maps of the SFR measured in each spatial bin within the time ranges specified in the upper left corner of each panel. At the youngest ages (upper left), we have grouped the time bins together into bins of roughly 25 Myr to create smoother maps for our analysis. Each analysis region is represented as a pixel, with the color of the pixel corresponding to the best-fit SFR during the time span specified. Regions plotted in black have a measured SFR of below 10 −4 M yr −1 , comparable to the typical uncertainty in low SFR regions. There are clear time-dependent changes in the pattern of SF, which we discuss in Section 5.
We note that towards the older bins in Figure 6, the central region of the galaxy does not have detectable star formation, suggesting that the significant star formation in the galaxy center may have started only ∼100 Myr ago. Our data in this region are sensitive to these stellar ages, so this timescale for the star formation at the galaxy center may be robust. However, this region also shows significant dust content (see Figure 8), with high values of dA V , suggesting that we may not be seeing the full characteristics of the young populations in this region.
We also present the integrated SFH over the last ∼630 Myr in Figure 7. To generate this figure, we summed the total stellar mass formed in each sub-region over all time bins back to 631 Myr. We find the integrated SFH shows clear spiral structure with a prominent two arm spiral that looks quite different from the more flocculent structure visible in an optical image. In addition there is a noticeable bar structure. We return to these points in detail in Section 5.

SFR Over Different Timescales
It is instructive to calculate the integrated SFR over various timescales (0−10, 0−100, 0−500, and 0−630 Myr) for comparison with different extragalactic SFR indicators, which we include in Table 5. We also include the SFR over these timescales for each individual region in Table 3.
The average SFR within the PHATTER footprint within the past 100 Myr is useful to compare with SFR measurements involving broadband FUV flux (see Kennicutt & Evans 2012). When restricted to the most recent 100 Myr, the average SFR within the PHATTER survey is 0.32±0.02 M yr −1 or 8.4±0.5×10 −3 M yr −1 kpc −2 . We plot the spatially-resolved SFR over the past 100 Myr in Figure 11, which we show alongside a FUV image of M33 ).

Recovered Extinction
We measure relatively modest foreground A V values across the PHATTER survey footprint with values between 0.0 and 0.5 mag. We measure a wider range in differential extinction, finding dA V values ranging from 0.0 to 2.5 magnitudes, with larger values of dA V in the center and concentrated in spiral arms. When we plot the dA V values for regions across the PHATTER footprint, we find that the morphology follows the galaxy's spiral structure, which suggests that we recover the intrinsic dust extinction in our fits. In Figure 8 we present a map showing the best-fit dA V value for each analysis region. We include the best-fit A V and dA V value for each region in Table 3.
For all regions with high measured dust extinction (A V + dA V > 2.3), we looked at each region's CMD and measured SFH by eye to confirm that the high dust extinction wasn't a result of a failure of our fitting code. There was only one region for which the high A V + dA V value was the result of a catastrophic failure of the code, region 1595 which contained a shredded foreground star. We physically masked the foreground star and its diffraction spikes out of the region's photometry and re-calculated its SFH. We present the corrected SFH for region 1595 in all results.
Even though our measured A V and dA V values seem to accurately recover the known foreground extinction to M33 and the expected dust morphology, we still expect to miss some star formation in M33 due to embedded star formation. The geometry of the dust and young stars -whether the stars are in front of or behind the dust from our vantage point -will affect our measurement of dA V and the star formation rate. We expect to miss some star formation in regions where the young stars are hidden behind too much dust. As shown in Figure 2, at the largest dA V values, some of the young stars are reddened all the way into the exclusion region of the CMD, meaning that our measurements are not sensitive to levels of dust reddening beyond ∼ A V + dA V 2.4. We expect that our measured 100 Myr SFR for a given region should correlate with the infrared emission from dust warmed by young stars (24 µm). To estimate the fraction of star formation that we miss due to embedded star formation, we look for regions where the 24 µm emission suggests that there is significant star formation, but our measured 100 Myr SFR is lower than expected. In Figure 9 we plot our measured 100 Myr SFR for each analysis region against the 24 µm flux in that region from Spitzer (Tabatabaei et al. 2007). We re-binned the Spitzer image to match the pixel scale of our SFH maps for comparison and plot relative units of flux within the re-binned image. The point for each region is color coded by its best-fit A V + dA V . We include the line of highest correlation as a black dotted line and include the 1σ range in gray.
The clear correlation between our rates and the warm dust emission suggests that we do capture most of the star formation associated with dusty regions. However, we expect to have underestimated star formation in the regions that lie far above the line of highest correlation in the left panel of Figure 9, as these regions have lower 100 Myr SFR than their 24 µm flux predicts, suggesting we have missed star formation embedded in dust.
To quantitatively look for potential embedded star formation, we visually inspected optical and 24 µm images of regions that lie above the 1σ range in Figure 9. We generated optical finder images using the PHATTER mosaic images in the following bands: F160W (red), F814W (green), and F475W (blue). We generate IR finder images using the full-resolution Spitzer 24 µm images of M33 (Tabatabaei et al. 2007) 9. In Figure 9 we present optical and IR finders for two different regions: one where the 24 µm flux suggests we may be missing embedded star formation (Region 978) and one where the 24 µm flux and 100 Myr SFR agree (Region 1825). In the optical image of Region 978, we see blue light from active star formation and a dust cloud that is most visible just below the bright blue star forming region at the top of the image. In the IR image, the white square outlines the area shown in the optical image. The 24 µm flux and 100 Myr SFR for this region are plotted on the left panel, outlined in a black star. We see that the 100 Myr SFR is lower for this region than we would expect given its bright 24 µm flux, suggesting that we may be missing some star formation embedded within the dust. The 100 Myr SFR measured for region 1825, on the other hand, matches up with the region's 24 µm flux (outlined with a black diamond on the left panel), suggesting that even though warm dust is present in this region, we are not missing significant star formation.
To quantitatively estimate the fraction of star formation missed, we can look at the outliers in the left panel of Figure 9. If we assume that at least some of the outliers are indicating missed star formation (for an example of this, at more than 2 sigma above the line), we can calculate the amount the SFR in each of these outlying regions would need to be increased to move them onto the line of highest correlation. We identify 21 regions that lie above 2σ and estimate that in these regions, we could be missing ∼7% of the total star formation within the last 100 Myr due to embedded star formation.
The results presented in this section describe the outcome of our A V and dA V fits with the Padova models. We also fit for extinction with the MIST models. When we compare the measured A V and dA V , we find that the two model sets agree on recovered extinction tile to tile with a mean difference in the A V measurements of 0.03±0.07 mag. The mean tile to tile difference in dA V is 0.01±0.3 mag, suggesting that there is not an offset in the recovered dA V values, but that there is scatter. There is not a spatial trend across the galaxy in the difference between the extinction measurements with the Padova and MIST models.

Recovered Metallicity
Our recent star formation histories should not be strongly sensitive to metallicity because the color of the upper main sequence CMD feature that our SFHs rely on is not strongly impacted by metallicity (e.g., Bressan et al. 2012). However, there is some sensitivity from the colors of more evolved Helium burning stars, which could potentially yield a signal when combined over many measurements. We did allow some freedom in the metallicity to be fitted between -2.3<[M/H]<0.1 for the Padova models and between -2.0<[M/H]<0.5 for the MIST models. One way to determine if our SFH results have any metallicity sensitivity is to compare the slope of the best-fitting mean metallicities with radius to the slope of the known M33 gas-phase metallicity gradient, such as that measured with observations of H II regions and planetary nebulae (e.g., Magrini et al. 2010).
To make this comparison we calculated the mean metallicity recovered in our SFH fits radially across the PHATTER survey footprint. We binned the SFH regions into annuli with a width of 1 kpc from the D25 ellipse center of M33 ) and measured the mean metallicity in each radial bin. Due to the large uncertainties on our metallicity measurements, we find that the slope of the metallicity gradient recovered with both the Padova and MIST models is consistent with both a flat gradient and the slope of the gas-phase metallicity gradient from Magrini et al. (2010), within 1σ errors. Additionally, we do not see a correlation (either positive or negative) between our measured metallicity and A V or our measured metallicity and dA V when we compare tile to tile, suggesting that we are not seeing metallicity as an artifact of the metallicity-extinction degeneracy.

DISCUSSION
The SFH maps reveal the time evolution of several interesting structures in the disk of M33, among these are the spiral arms and a central bar. We discuss both of these examples in detail below.

Evolution of Spiral Structure
The time resolved structure in the SFR maps in Figure  6 show clear spiral structure. However, this structure appears qualitatively different at early and late times. In the earliest time bins (top row), we observe a flocculent spiral structure between the present day and ∼ 79 Myr ago. Starting in the fourth time bin (79−100 Myr ago) however, we begin to see the emergence of a strong two armed spiral structure, which is evident in all of the time bins >100 Myr. This change in morphology can be difficult to see at the full time resolution of Figure  6. To better visualize the evolution in spiral structure and reduce the impact of random small scale fluctuations in SFR measurements, we combined the SFR measurements for all bins younger than 79 Myr (which are somewhat noisy due to the stochastic nature of SF over short time intervals) and all bins older than 79 Myr. We present these time-integrated SFH maps in Figure 10, showing the pronounced difference between the structures formed by younger and older stars.
M33 has typically been characterized as a flocculent spiral galaxy with a number of weaker spiral arms, rather than the strong arms seen in grand design spirals (e.g., Humphreys & Sandage 1980;Dobbs et al. 2018). The spiral structure in the left panel of Figure 10 (0−79 Myr ago) shows multiple spiral arms, in agreement with the morphology that M33 shows in optical images. This correspondence suggests that the optical morphology is strongly shaped by the very low mass-to-light (M/L) ratio populations of the youngest stars. However, it important to note that this multi-arm spiral does seem to include the two spiral arms seen in the two-arm spiral structure shown by the morphology of the older stellar populations (older than 79 Myr, right panel of Figure  10).
M33's lack of optically visible strong spiral arms has led to some debate in the literature about whether the spiral structure of M33 is a result of tidal interactions with M31 (e.g. Semczuk et al. 2018) or if this spiral structure arose without interaction (van der Marel et al. 2019), as a result of gravitational instabilities in the stars and gas in M33 (Dobbs et al. 2018). Near infrared (NIR) observations in the J, H, K bands and with 2MASS have previously revealed a stronger two-arm spiral structure within the inner disk of M33 (Regan & Vogel 1994;Jarrett et al. 2003).
The two-arm structure seen in the NIR agrees with our map of the SFH of M33 between 79 and 631 Myr ago. On the surface, this agreement may be surprising given the differing timescales of the two tracers. The NIR is typically presumed to be dominated by older red giant branch (RGB) and asymptotic giant branch (AGB) stars, far older than 500 Myr. However, there are two possible explanations. First, the NIR is in fact sensitive to red core Helium burning stars, which appear for ages 300 Myr (Regan & Vogel 1994;Melbourne et al. 2012). Second, a strong long-lived spiral pattern in the stellar mass (as traced by the RGB stars that dominate the disk mass) can act as a gravitational perturbation that helps confine the younger stars, following global density wave theory (e.g., Lindblad 1960;Lin & Shu 1964). Thus, while the gravitational-instability driven flocculent structure in the gas is imposed on stars at birth, this structure dissipates over ∼50 Myr as stars migrate from their birthplace and older stars are swept into the larger tidally-driven gravitational perturbation (e.g., Dobbs & Pringle 2010).
A more detailed study of the structure of M33 with resolved stellar populations from the new PHATTER dataset is forthcoming, including a study of the age gradient across the spiral arms (A. Smercina et al., in preparation).

Detection and Formation of M33 Bar
There has been some debate in the literature about whether M33 hosts a central bar. A weak central bar has been observed via photometric measurements at infrared (Regan & Vogel 1994;Block et al. 2004;Hernández-López et al. 2009) and blue wavelengths (Elmegreen et al. 1992), and via measurements of gas dynamics (Corbelli & Walterbos 2007). However, Regan & Vogel (1994) debated whether the weak central bar is actually an inner extension of the two spiral arms. Subsequent work studying gas dynamics in the center of M33 suggests that a small bar could explain discrepancies between models and the observed gas velocities in the inner (r<2kpc) disk (Kam et al. 2015). Theoretically, Sellwood et al. (2019) and Semczuk et al. (2018) have both presented simulations that favor the formation of a bar in the inner 3 kpc of M33 within about 1 Gyr. For reference, the PHATTER survey extends ∼3.5 kpc from the center of M33, and ∼2 kpc along the semi-major axis.
While the young populations near the center of M33 appear to be confined to a relatively small and circular region, we observe a bar structure in the stellar distribution of stars older than roughly 79 Myr. This bar is also clearly visible in the integrated SFH map of star formation back to 630 Myr (Figure 7). The maps of the integrated SFR in Figure 10 show no trace of the bar in populations younger than 79 Myr. The panel on the right, showing the SFR between 79 and 631 Myr ago, shows a clear bar extending out to ∼1 kpc. This suggests that the gas (which sets the spatial distribution of the youngest stars), does not trace the structure of the stellar bar, raising an interesting dynamical question about why the gas is able to resist formation of a bar, while the older stellar population is not.
More detailed analysis to fully characterize observations of the bar in M33 with the PHATTER data will be included in future work on the galaxy structure (A. Smercina et al., in preparation), but its appearance in the star formation history maps potentially gives us a timescale for stars to be perturbed into density waves and bar like orbits.

Comparison of SFH results with other broadband measurements
We find that the morphology we observe in our 0−100 Myr SFH maps follows the morphology visible at FUV wavelengths. In Figure 11, we present a FUV image of M33 ) alongside a map of our SFR measurements in each spatial region over the past 100 Myr. Regions with high SFR within the past 100 Myr spatially coincide with high flux areas in the FUV image.
The PHATTER survey only covers the inner disk of M33, so to quantitatively compare our SFR measurements with published values, we have to correct our measured SFR for the star formation missed outside our survey area. For consistency with published work, we correct to an elliptical aperture within the D25 isophote, using the FUV D25 ellipse from Gil de . We then measure the total FUV flux within this D25 ellipse and determine the fraction of the flux contained within the PHATTER survey footprint. We find that the fraction of total flux within our survey footprint is ∼45%, meaning we need to multiply our total SFR measurement by a factor of ∼2.2 to compare with global FUV SFR measurements. When we scale up our SFR by this factor, we get a total SFR of ∼0.74 M yr −1 for M33 within the last 100 Myr.
This value can be compared to the analysis of Verley et al. (2009), who measured the global SFR of M33 using a combination of Hα and FUV flux measurements, using IR luminosities to correct for dust extinction. Their extinction corrected value for the global SFR in M33 over the last 100 Myr is 0.45±0.10 M yr −1 . This measurement is lower than our scaled-up global SFR by a factor of about 1.6 and does not agree within errors, suggesting that the integrated light measurements may be missing some star formation that the resolved stars are able to capture. This underestimation agrees with findings by Lewis et al. (2017), who found that the 24 µm corrected FUV flux underestimated the SFR by a factor of 2.3−2.5 when compared to CMD-based SFHs in M31.
To check the correspondence between our SFHs and FUV emission more closely, we compared the SFH measured for some of our individual spatial regions that have high measured values of FUV and Hα flux. In Figure  12 we present the SFH for three regions: one with relatively FUV flux, one with relatively high Hα flux, and one with neither. We include a red-blue image where the red band is an Hα image of M33 from the Local Group Galaxy Survey by Massey et al. (2007)  The SFHs measured for these individual regions match our understanding of Hα and FUV flux as tracers of recent star formation. We note that most regions that had high Hα flux also had fairly significant FUV flux and visa versa. The region selected for high Hα flux shows a peak in star formation within the last 5 Myr, but also has consistently high star formation over the past 100 Myr. The region selected for high FUV but low Hα flux shows high star formation within the past 100 Myr, but does not show significant star formation within the past 5 Myr. The region that does not have significant FUV or Hα flux shows essentially no star formation within the past 100 Myr and very low rates of star formation between 100 and 630 Myr, out to the oldest age considered in our analysis. This excellent correspondence suggests there is much more to be learned from the detailed comparisons of our SFH measurements and the UV emission (and other integrated light) characteristics of regions in M33, including comparisons with both local extinction and age. We plan to complete such comparisons in future work.

CONCLUSIONS
We have measured the SFH of M33 within the opticalonly PHATTER survey footprint, which covers a deprojected area of ∼38 kpc 2 and extends to 3.5 kpc, 1−2 scale lengths, from the center of M33 . To measure the spatially resolved recent star formation history, we divided the survey area into 2005 analysis regions measuring 24 (roughly 100 pc) on a side and measured the SFH for each region using colormagnitude diagram fitting. We summarize our main findings below: 1. We measure the mean SFR over the last 100 Myr within the PHATTER survey area to be 0.32±0.02 M yr −1 . When we scale up our measurement to account for the full D25 area of M33, we measure a mean SFR over the last 100 Myr of ∼0.74 M yr −1 , which is higher than previous FUV 100 Myr SFR measurements.
2. We measure the current SFR (within the past 10 Myr) within the PHATTER survey area of M33 to be 0.20±0.03 M yr −1 .
3. We observe evolution in the spiral arm structure of M33 over time. The most recent SFH maps of M33 show a flocculent multi-armed structure, while the maps older than ∼ 79 Myr ago show a more distinct two-armed structure.
4. We observe a bar in the center of M33, which is observed in stellar populations older than ∼ 79 Myr. Younger stellar populations formed in more recent star formation episodes do not appear to follow the bar structure.
We thank Leo Girardi for his thoughtful comments on this paper. Support for this work was provided by NASA through grant #GO-14610 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. Margaret Lazzarini was supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102721 during a portion of this work. The difference between the best-fit SFH measurements using the Padova and MIST stellar models is plotted as a function of time. The dashed horizontal line marks a difference of 0. bottom: The black histogram represents the SFH when we use the Padova stellar models in our CMD fitting and the pink histogram shows the SFH when we use the MIST stellar models. The error bars represent the random uncertainties. The horizontal dashed lines mark the mean SFR value for each set of stellar models, the mean SFR values agree so closely, they are difficult to distinguish on the plot. The two SFH fits illustrate systematic uncertainties in our measurements. The models agree very well up to about 20 Myr ago. Differences at younger ages reflect differences between how the two model sets treat the most massive, young stars. For a full discussion of the differences between these two SFH measurements and systematic uncertainties in our SFH measurements, see Section 3.3. The SFR values and error bars plotted for the Padova stellar models are listed in Table 4. Due to the extent of our survey area, we expect to have measured roughly 45% of the SF in M33, which we determined by comparing to the fraction of total FUV flux within the PHATTER survey footprint.   Figure 6. SFH maps of M33 within the PHATTER footprint area over indicated timescales. These time bins do not reflect the full 0.1 log(year) resolution of our SFH fits, but we group the SFH up to 100 Myr ago into ∼25 Myr bins to more clearly show the structure of the spiral arms. The SFR values shown are the best-fit values from our fits. Each pixel represents a 100 pc by 100 pc region for which we derived the SFH using CMD-fitting. We note that the central pixels in the later time bins (older than 100 Myr) show 0 SFR, as noted in Section 4. The cause for this is not yet clear, but could be due to incompleteness.  Figure 7. The spatially resolved star formation rate of M33 between 0 and 631 Myr ago. This map is essentially a combination of the multiple SFH maps shown in Figure 6 across temporal bins. The SFR value in each bin was calculated by summing the stellar mass formed in that spatial bin over all time bins between 0 and 631 Myr ago and then dividing by the total elapsed time. Each pixel represents a 100 pc by 100 pc region for which we derived the SFH using CMD-fitting.  Figure 8. We present a map showing out best fit differential extinction (dAV ) value for each analysis region in the PHATTER survey area. We observe that the high extinction regions follow the spiral structure of M33, suggesting we adequately recovered dust extinction in our fits. We expect that our dAV maps may differ slightly from maps of dust emission due to the geometry of stars and dust in individual regions (i.e., if stars are in front of or behind dust clouds from our point of view). We describe our method for measuring dAV in Section 3.2.  (Tabatabaei et al. 2007). We plot the line of highest correlation with the black dotted line and show the 1σ range shaded in gray. The points are color coded based on the measured AV + dAV in that region. We expect to have under-estimated the star formation rate for regions that lie above the 1σ region, because those regions have lower 100 Myr SFR than expected based on their 24 µm flux, suggesting that we have missed embedded star formation. right, top row: We present an optical color image from the PHATTER survey and a 24 µm image from Spitzer for region 978 (outlined on the left panel with a black star), where we expect to have underestimated star formation due to embedded star formation. The PHATTER image measures 24 on a side and was created using the PHATTER survey mosaic images in the following bands: F160W (red), F814W (green), F475W (blue     (Massey et al. 2007). We highlight three regions and plot the best-fit SFH for those regions in the main panel. We selected a region with high FUV flux, a region with high Hα flux, and a region without high flux in either bandpass. The region without high flux in either FUV or Hα shows little to no star formation in the last 100 Myr and low star formation rates beyond 100 Myr. The region selected for high Hα flux shows significant star formation in the most recent time bin and high star formation rates back to about 100 Myr. The region selected for high FUV flux shows high star formation over the last 100 Myr but does not have measured star formation in the most recent time bin. These findings agree with our understanding of Hα emission as a tracer of very recent (within the past 5 Myr), or instantaneous star formation (e.g.; Byler et al. 2017) and FUV emission as a tracer of star formation within the last 100−200 Myr as it directly traces emission from young, massive stars (e.g.; Kennicutt & Evans 2012). As we would expect, a region without significant Hα or FUV emission does not show significant star formation within the past 100−200 Myr.  Dec. 4 SFR 6.6-6.7 SFR 6.7-6.8 SFR 6.8-6.9 SFR 6.9-7.0

a) b) c) d)
(1) (3) (8) (10) SFR 7.0-7.1 SFR 7.1-7.2 SFR 7.2-7.3 SFR 7.3-7.4 SFR 7.4-7.5 SFR 7.5-7.6 SFR 7.6-7.7 SFR 7.7-7.8 SFR 7.8-7.9 SFR 7.9-8.0 SFR 8.0-8.1 SFR 8.1-8.2  Note-We only include the SFH for the first 6 spatial regions in the print version of this paper. This table is available in its entirety in machinereadable form in the online version of this paper. See Table 2 for a description of each column. Note-The spatially-integrated SFH of the PHAT-TER survey in M33. The first column lists the start of each time bin in log(years) and the second column lists the end of each time bin in log(years). The third column lists the SFR for that time bin with uncertainties. The listed uncertainties incorporates both random and dust uncertainties. For more information on how these uncertainties were derived, see Section 3.3.

A. COMPARING MODEL CMDS WITH OBSERVED CMDS
In our fitting process, we generate model CMDs which can be compared against the observed CMDs to assess the quality of our fits. For each region being fit, we can then divide the CMD into bins and compare the number of observed stars in each bin to the number of model CMD stars in each bin. We can then find the difference between these two values and calculate the significance of this difference. The significance is calculated using a Poisson maximum likelihood statistic, based on the Poisson equivalent of χ 2 (Dolphin 2002). To demonstrate how well both the Padova and MIST models perform in creating model CMDs that fit the data, we have included figures (13, 14, and 15)  shows the difference between the observed and Padova model CMD and the middle right panel shows the significance of the difference, which is described above. The lower left panel shows the difference between the observed and MIST model CMD and the lower right panel shows the significance of the difference. We can analyze the significance of the difference between the model and observed CMDs to look for 1. differences between how well the MIST and Padova model CMDs fit the data and 2. trends that could suggest systematic bias in our SFH measurements. In Figures 13,  14, and 15 there is not a clear difference in the significance of the difference between the model and observed CMDs for the MIST and Padova models, which we can see when comparing the middle right and lower right panels of these figures. We also do not see a strong trend in the significance of the difference between the model and observed CMDs in the areas of the CMD being fit to recover the SFH. Specifically, the difference between the model and observed CMDs looks fairly even across the CMD, rather than demonstrating that our models under or over produce stars in different parts of the CMD that are outside of the shaded region excluded in our fits. The significance of the difference between the model CMDs and real CMDs is used to calculate the goodness of the fit, which is based on the Poisson equivalent of χ 2 . To get a holistic view of how well the Padova and MIST models represent the whole data set, we compared the fit values for each of the 2006 regions in our SFH maps. We found that the distribution of fit values across all regions was very similar for the Padova and MIST models. The Padova mean fit value across all regions agrees with the MIST mean fit value across all regions within <1%. This suggests that while there are some differences at young ages (<20 Myr) between the Padova and MIST models due to theoretical uncertainties, both sets of stellar models are similarly successful at matching the data.

B. SFH RECOVERY TESTS
To assess how well our CMD fitting performed, we generated mock CMDs using known input SFH and extinction values, A V and dA V values. We can then compare the SFH and extinction that are measured in the same way we measure these values for real data. We performed two sets of tests. First, we used realistic SFHs to create mock CMDs. We also ran tests where we used simplistic SFHs with a burst of star formation in a single bin to assess how well the bin and strength of these bursts were recovered.

B.1. Realistic SFH Tests
For a given test, we chose a region from our SFH maps and used the SFH that was measured for that region as the input to generate a mock CMD. We used that region's best-fit A V and dA V when creating the mock CMD. We then ran this mock CMD through the same pipeline we used to measure our SFHs from the data. We performed this test for regions of varying stellar density and extinction. We chose regions that had combinations of high/medium/low stellar density and high/medium/low extinction. We did not perform these tests for high stellar density/low extinction or low stellar density/high extinction because we did not find any regions in our data that showed this combination of stellar density and extinction, thus these tests would not be representative of our data. We performed these tests with both the Padova and MIST models, using the same model to create the mock CMD and fit for the SFH. We compared the input versus output A V and dA V values for our test regions. We calculated the root mean square error, which is The difference between the observed and MIST model CMD (left) and the residual significance (right). We find no systematic residuals in the area of the CMD included in our fits, indicating that both sets of stellar models are a good fit to the data. The difference between the observed and Padova model CMD (left) and the residual significance (right), which is the observed CMD minus the modeled CMD, weighted by the variance. Bottom: The difference between the observed and MIST model CMD (left) and the residual significance (right). We find no systematic residuals in the area of the CMD included in our fits, indicating that both sets of stellar models are a good fit to the data. MIST Models Figure 16. Comparison between the input and output AV and dAV values tests where we generated mock CMDs and recovered their extinction and SFH. We created mock CMDs with the input AV and dAV values and then measured these values using the methods described in Section 3.2. We used the SFH and extinction values for regions from our SFH maps that had varying stellar density and extinction. Filled points indicate AV values and open points indicate dAV values. Dark purple points are high extinction regions, pink points are medium extinction regions, and orange points are low extinction regions. Circles are high stellar density regions, stars are medium stellar density regions, and squares are low stellar density regions.
defined as the square root of the sum of the squares of the differences between the input and output values, scaled by the total number of data points. We found that for the Padova models, RM SE Av = 0.071 and RM SE dAv = 0.25. For the MIST models, RM SE Av = 0.093 and RM SE dAv = 0.13. In Figure 16 we plot the input versus output A V and dA V values for our tests with both model sets. Filled points indicate A V values and open points indicate dA V values. Dark purple points are high extinction regions, pink points are medium extinction regions, and orange points are low extinction regions. Circles are high stellar density regions, stars are medium stellar density regions, and squares are low stellar density regions. Both model sets do a good job of recovering A V and dA V values, with the highest extinction regions presenting the most challenge in accurately recovering the dA V value. We also compare the input and output SFH for these test regions. There are several ways to make this comparison including calculating the median age of the region and comparing the cumulative fraction of stellar mass formed over the time range analyzed in our fits (0-630 Myr ago). We compared the median age of the SFH used to create the mock CMDs against the median age of the SFH measured with these mock CMDs for both the Padova and MIST models. In this analysis, we use the SFH within the last 630 Myr, the age back to which our fits are reliable. The median age is the age at which 50% of the total stellar mass formed, and the lower and upper errors come from the ages at which 16 and 84% of the stellar mass formed. The input and output median ages agree within errors for both model sets, however the Padova models do a slightly better job recovering the median ages for regions with median ages older than about 300 Myr. We can also compare the cumulative stellar mass fraction formed to compare the input and output SFH over the last 630 Myr. We show these cumulative mass fractions for the input SFH and the output SFH retrieved with the Padova and MIST models in Figures 17 and 18 for regions of varying stellar density and extinction. In Figure 17 all regions have roughly the median value of A V and dA V in our dataset, but vary in their stellar densities. In Figure 18, all regions have roughly the median stellar density of regions in our dataset, but have varying levels of extinction. The errorbars on the output cumulative stellar mass fraction in both figures were generated by randomly sampling within the errors on the SFR in each time bin 1000 times. low stellar density medium extinction Figure 17. Comparison of the cumulative fraction of stellar mass formed over time from tests where we recovered the SFH from mock CMDs generated with realistic SFHs and varying stellar density. The orange line represents the SFH that was used to generate mock CMDs for these regions. The dark purple line represents the SFH recovered from the mock CMD that was fit with the Padova models and the pink line represents the SFH recovered from the mock CMD that was fit with the MIST models. This figure demonstrates the effect of stellar density on our measurements. For the high stellar density region with medium extinction, our errors are fairly small because of the large number of stars in each region on the CMD. The errors are larger, by contrast, in the low stellar density region where the errors are dominated by statistical errors stemming from the low number of stars on the CMD. medium stellar density low extinction Figure 18. Comparison of the cumulative fraction of stellar mass formed over time from tests where we recovered the SFH from mock CMDs generated with realistic SFHs and varying levels of extinction. The orange line represents the SFH that was used to generate mock CMDs for these regions. The dark purple line represents the SFH recovered from the mock CMD that was fit with the Padova models and the pink line represents the SFH recovered from the mock CMD that was fit with the MIST models. This figure demonstrates the effect of extinction on our measurements. The SFH is recovered quite well by both models in all three regions, with smaller errors in the region with low extinction.
within the errors. In summary, we ran tests where we created mock CMDs using a realistic SFH and varying levels of stellar density and extinction. We found that both models recovered the level of extinction quite well, with the MIST models performing slightly better at higher levels of differential extinction. When comparing the median ages of the input SFH versus measured SFH, the Padova models did a slightly better job recovering the median ages towards the older half of our fits. We also compared the cumulative stellar mass fraction for the input SFH and SFHs measured from our mock CMDs. We found that there was not a discernible difference in how well the MIST and Padova models recovered the cumulative stellar mass fraction. We find the largest difference between the input and output SFH for regions with low stellar density, reinforcing the notion that it is better to have a large number of stars on the CMD when recovering the SFH. Note-Summary of tests where a mock CMD was generated with a burst SFH, which only had SFR in one time bin. Column 1 lists the name of the test. Column 2 lists the time bin of the burst in log(yr) and column 3 lists the start time of the burst in Myr. Column 4 indicates the SFR during the burst. Columns 5−7 indicate the fraction of the 100 iterations of a given burst tests during which the burst was recovered in the exact input time bin, within one bin of the input time bin, and within two bins of the input time bin. Column 8 lists the burst strength fractional difference, which is defined as (output burst SFR − input burst SFR)/input burst SFR.