Coevolution of Stars and Gas: Using an Analysis of Synthetic Observations to Investigate the Star–Gas Correlation in STARFORGE

We explore the relation between stellar surface density and gas surface density (the star–gas, or S-G, correlation) in a 20,000 M ⊙ simulation from the STAR FORmation in Gaseous Environments (starforge) project. We create synthetic observations based on the Spitzer and Herschel telescopes by modeling contamination by active galactic nuclei, smoothing based on angular resolution, cropping the field of view, and removing close neighbors and low-mass sources. We extract S-G properties such as the dense gas-mass fraction, the Class II:I ratio, and the S-G correlation (ΣYSO/Σgas) from the simulation and compare them to observations of giant molecular clouds, young clusters, and star-forming regions, as well as to analytical models. We find that the simulation reproduces trends in the counts of young stellar objects and the median slope of the S-G correlation. This implies that the S-G correlation is not simply the result of observational biases, but is in fact a real effect. However, other statistics, such as the Class II:I ratio and dense gas-mass fraction, do not always match observed equivalents in nearby clouds. This motivates further observations covering the full simulation age range and more realistic modeling of cloud formation.


INTRODUCTION
The majority of stars form in associations or groups within giant molecular clouds (GMCs, Lada et al. 1991;Krumholz et al. 2019;Cheng et al. 2022), which can vary greatly in size, from ∼10 to thousands of stars (Porras et al. 2004).Feedback from embedded clusters often quickly disperses the natal clump or even the entire GMC (Lada 2005;Krause et al. 2020).Therefore, the relationship between gas and young stellar object (YSO) density provides important clues about the star formation process and cloud evolution.Schmidt (1959) was one of the first to present an analytical model of the relationship between star formation rate (SFR), and thus stellar mass, and gas density.That work suggested that SFR and gas density follow a power law relationship.
This correlation was examined over the next several decades by a number of authors (e.g.Sanduleak 1969;Hartwick 1971).However, it was not until improved observational capabilities and analysis techniques in the 1980s and 1990s (e.g.Kennicutt 1989Kennicutt , 1998) that strong evidence was found for its viability.This work motivated an analogous relation known as the Kennicutt-Schmidt (KS) law that applies to line-of-sight surface densities of gas and the star formation rate per area: Henceforth, we refer to this relation as the Star-Gas or S-G correlation.This relationship has since been wellcharacterized as a power-law with an index of N∼1.4 as applied to galaxy-scale star formation (see Kennicutt & Evans (2012) for a detailed review).
At smaller scales within individual galaxies, there is also evidence for the presence of an S-G correlation.For example, Bigiel et al. (2008) used HI, CO, 24 µm, and UV data to examine the S-G correlation at 750 pc resolution in 18 nearby spiral and dwarf galaxies.Many regions showed a strong power-law relation, although the power-law index varied from 1.1 to 2.7 based on position.They also observed that the star formation efficiency (SFE) decreased with galactic radius, which they argued implies a connection between environment and the S-G correlation.
However, the methods used to measure the SFR on ≳ kpc scales, such as Hα, far-UV, and 24 µm emission, become less effective at smaller spacial scales.The results of Liu et al. (2011), as well as modeling by Calzetti et al. (2012) show that this kind of analysis breaks down with shrinking sample area because star formation is not well-sampled statistically.Gutermuth et al. (2011) (G11 hereafter) demonstrated that the SFR calculated from far-IR luminosity (L FIR , e.g., Heiderman et al. 2010) underestimates the SFR calculated from counts of YSOs in nearby young clusters by up to an order of magnitude.This is because measurements based on far-IR luminosity assume a well-sampled stellar initial mass function (IMF) and reliable sampling of the GMC mass function to fully sample the lifetime of high-mass stars.However, in order to satisfy these assumptions, measurements must be integrated over physical scales ≳ 1 kpc (Calzetti et al. 2012).
To avoid the smoothing inherent to measurements of star formation relations in other galaxies, some recent studies instead focus on individual star-forming regions in the local Milky Way, where it is possible to identify and count individual forming stars with high completeness.Since YSOs provide a direct measurement of the SFR, a simple estimate of the total mass converted to stars per time is given by where m YSO is the average mass of a YSO, n YSO is the number of YSOs, and t avg is the characteristic timescale for the YSO evolutionary stage or stages considered.
By utilizing YSO censuses from Spitzer, G11 and Pokhrel et al. (2020) (P20 hereafter) found and measured an intracloud S-G correlation with an index of N ≈ 2 in several nearby GMCs.While initial measurements varied widely (N = 1.5 -4) (G11), P20 reduced intrinsic scatter in the measurements by adopting a uniform YSO extraction from the Spitzer Extended Solar Neighbor Archive (SESNA), utilizing more robust Herschel -based GMC gas column density maps, and by specifically using YSOs in the early stages of star formation.This led to N = 1.8 -2.3 in 12 nearby clouds with gas masses varying over three orders of magnitude.Also, the scaling factor in the S-G correlation varies between clouds (Lada et al. 2013, G11, P20), but the scatter in the scaling factor is reduced significantly when it is normalized by the gas freefall time (Pokhrel et al. 2021).This implies that the SFE per freefall time has limited variation, which may indicate that local processes (e.g., protostellar outflows and stellar winds) govern and regulate star formation (Guszejnov et al. 2021;Pokhrel et al. 2021;Hu et al. 2022).
In order to gain a better understanding of how local processes impact star formation, it is useful to turn to theoretical models and numerical simulations.However, observed S-G correlations have only recently started to become incorporated as constraints for models of starforming molecular gas.P20 used simulations by Qian et al. (2015) that used the ORION adaptive mesh refinement code (Truelove et al. 1998;Klein 1999) to create synthetic observations similar to observations taken by Herschel.That work reproduced similar S-G correlations for 12 nearby GMCs using hydrodynamic turbulent simulations and an analytical model of thermal fragmentation.While the simulation produced an S-G correlation that is very similar to observations, it did not include magnetic fields or kinematic feedback.In this work, we analyze a 20, 000 M ⊙ run of the STAR FORmation in Gaseous Environments (starforge) project, the first massive GMC magnetohydrodynamics simulation to resolve individual stars while including multiband radiation, stellar winds, protostellar outflows, and supernovae (Grudić et al. 2021(Grudić et al. , 2022, etc.), etc.).
In order to most effectively compare the starforge simulation to observations, we construct synthetic observations according to the data used in P20, taking into account the known specifications and limitations of Spitzer and Herschel data.In Section 2, we describe the specifics of the simulation snapshots and our methods for creating synthetic observations.In Section 3, we present results from our investigation into various star-gas properties in the simulation and compare to observations.Discussion is provided in Section 4, and a summary and conclusions are given in Section 5.

starforge Simulations
The starforge framework is built on the gizmo meshless finite mass magnetohydrodynamics code (Hopkins 2015).The framework includes a variety of modifications that enable the modeling of individual forming stars and the interactions that occur with their cloud environment.In this work we analyze the starforge simulation presented in Grudić et al. (2022).We briefly summarize the simulation properties here and refer the reader to Grudić et al. (2021) for a detailed description of the starforge numerical methods.
The simulation follows the evolution of a 20,000 M ⊙ cloud with initial radius of 10 pc.The cloud turbulence is initialized so that the cloud is virialized with α ≡ 5σ 2 3D R cloud /(3GM cloud ) = 2, where σ 3D is the gas velocity dispersion.The initial magnetic field is uniform in the ẑ direction and corresponds to a massto-flux ratio relative to the critical value for stability µ ≡ 0.4 E grav /E mag = 4.2, where E grav and E mag are the total gravitational and magnetic energies, respectively.
The calculation follows the gas thermodynamics selfconsistently, including treatment of line cooling, cosmicray heating, dust cooling and heating, photoelectric heating, hydrogen photoionization, and collisional excitation of both hydrogen and helium.The evolution of the dust temperature is coupled to the radiative transfer step.gizmo's radiation transfer module follows five bands, which cover the frequencies corresponding to ionizating radiation, FUV, NUV, optical-NIR, and FIR (Hopkins & Grudić 2019;Hopkins et al. 2020).
Once gas satisfies multiple criteria intended to identify centers of unstable collapse, Lagrangian sink particles are inserted, which occurs at densities of ρ max ∼ 10 −14 g cm −3 .
The cell mass resolution is dm = 10 −3 M ⊙ , which allows the calculation to resolve the stellar mass spectrum down to ∼ 0.1 M ⊙ .The sink particles, henceforth referred to as stars, follow a sub-grid model for protostellar evolution and radiative feedback as described in Offner et al. (2009).The particles are also coupled to models describing protostellar outflow launching, stellar winds, and supernovae (Cunningham et al. 2011;Guszejnov et al. 2021;Grudić et al. 2021).The calculation continues until stellar feedback disperses the natal cloud and star formation concludes, which happens at ∼ 9 Myr.
The simulation has a final SFE of 8% that agrees with statistical models of nearby galaxies.Protostellar jets dominate feedback for most of the simulation and are important for regulating the IMF, but they cannot wholly disrupt the cloud.Eventually, radiation and winds from massive stars create bubbles that expand and disrupt the cloud, drastically reducing SF.By following GMC evolution, Grudić et al. (2022) measure a relatively unambiguous IMF.It resembles the Chabrier IMF with a high-mass slope of α = −2±0.1.The IMF is much more realistic than previous simulations without full feedback.Feedback from radiation/winds of massive stars limits the maximum observed mass to 55 solar masses, moderating the high-mass tail of the IMF.The integrated luminosity and ionizing photon rate are also very close to an equal-mass cluster with a canonical IMF.A more detailed study of the impact of various feedback processes and cloud initial conditions on the IMF is presented in Guszejnov et al. (2022).Grudić et al. (2022) also note the importance of directly comparing observations and simulations via synthetic observations, as we aim to do in this work.
To construct the stellar surface density, we require a minimum of 11 YSOs.The first snapshot with at least this number of sources is at 1.47 Myr.Altogether our analysis uses 16 snapshots, spaced 0.49 Myr apart, which span 1.47 to 8.80 Myr.

Constructing Synthetic Observations
For our analysis to better mirror that of P20, we create synthetic observations by including various considerations to bring our data closer to that which might have been observed by Spitzer and Herschel.We refer to analysis done with minimal adjustments, i.e., only 2D projection, age-to-class conversion, and 0.01 M ⊙ mass cutoff (see below) as the "fiducial analysis", while analysis with further considerations are collectively referred to as "synthetic observations."The fiducial (minimallyadjusted) case allows us to examine how well the simulation can reproduce various statistics and identify where observational biases may affect the agreement.In order to create these synthetic observations, we extract or compute the (line-of-sight-projected when applicable) molecular number density of H 2 and the masses, coordinates, ages, and particle indexes of the sink particles, which represent YSOs.

YSOs
YSOs fall into distinct groups based on their observed properties.Historically, these have been binned into representative classes (Lada 1987;Shu et al. 1987;Greene et al. 1994;Robitaille et al. 2007;Dunham et al. 2015), e.g., Class I, Class II, and Class III.1 .Note that class does not have a direct mapping to source age, but it is often used as a proxy for evolutionary stage.YSOs in each class differ in the shape of their spectral energy distribution (SED), which depends on the characteristics of the circumstellar material around the YSO.Class Is are usually deeply embedded in cold, dense, and dusty gaseous envelopes, Class IIs have classical protoplanetary disks, and Class IIIs have mostly lost their disks (or the visible disk material has substantially coalesced into larger planetesimals that are generally invisible in the infrared).
For the first step of our analysis, we map each of the starforge stars to an observational class.Ideally, the stellar age would be employed to directly map each source to the appropriate spectral class.However, the average age and lifetime of each class is uncertain, since the individual classes are not completely distinct and the boundaries between them are somewhat arbitrary.Class lifetimes are inferred observationally using the relative number of sources in each class and by assuming a typical disk lifetime (e.g.Dunham et al. 2014).Consequently, a full self-consistent class assignment requires constructing synthetic observations using radiative transfer to model the SEDs.Instead, we assign each star to a class based on its age (the time elapsed since the sink particle forms in the simulation) and adopt a statistical approach rather than an exact mapping.
We model the transitions from Class I → Class II and Class II → Class III as exponential decays, adapting the models and half-lives of the transitions from Kristensen & Dunham (2018) and Mamajek (2009) to represent the age-to-class conversion.Using these half-lives, we calculate two numbers corresponding to each source, f a and f b , which corresponds to the statistical weighting given to each source for transitions (a) (Class I → II) and (b) (Class II → III): where t age is the age of the YSO and t 1/2;a and t 1/2;b are the half-lives of the Class I → Class II transition (0.22 Myr) and the Class II → Class III transition (1.7 Myr), respectively.Then, we generate two random numbers (r a and r b ) for each source using consistent seeds and the persistent source index from starforge, so that each YSO has the same r a and r b for the entire run.If r a < f a , then the YSO is assigned to Class I.If not, we check whether r b < f b .If so, the YSO is assigned to Class II.If not, it is assigned to Class III.
By fixing r a and r b for each source, we ensure that the sources progress forward through the classes (I to II to III) as they age in the simulation.However, in actual observations, a YSO's trajectory may not be so linear.For example, Dunham et al. (2010) used models of accreting sources to show that YSOs undergoing episodic mass accretion may transition to an earlier Class.The notion that older sources can populate the earlier classes is also supported by the work of Hernández et al. (2007a), who observed what appear to be older, "evolved" disks.Another problematic assumption is that the Class lifetimes are the same for every environment, which is unlikely since protostars in areas of high YSO density tend to have greater luminosity (Kryukova et al. 2014;Cheng et al. 2022).Despite the approximate nature of our model for Class assignment, we find that it reproduces the expected YSO distributions well.Whereas, assuming an exact one-to-one mapping between age and Class leads to sharp transitions that do not match observations as closely.
Next, in order to model source confusion present in Spitzer observations, we inject Active Galactic Nuclei (AGN) contaminants.In Spitzer observations, background AGN can appear as YSOs of Class I and II with roughly equal probability (Gutermuth et al. 2008(Gutermuth et al. , 2009)).To simulate this effect, we randomly place N Class Is and IIs within the dataset, where N was determined to be ∼9 per square degree (P20).This has the immediate effect of introducing many sources with low spatial density.This is especially significant for the synthetic clouds at closer distances due to the commensurately larger angular size of the cloud (see Figure 1, where it is clear that AGN dominate over YSOs in low gas density regions).We then correct for these contaminants following the method used by G11: we adopt a threshold of log Σ gas > 1.3 M ⊙ pc −2 for points on the S-G plot (see Section 3.4).We adopt the same distribution of AGN contaminants for all snapshots to ensure that the AGN stay the same (i.e.same position and Class).
We also model instrumental detection limits to account for undetectable low-luminosity sources.To replicate this in the synthetic observations, we implement a simple mass cutoff, where we remove sources below 0.08 M ⊙ (200 and 400 pc distance) or 0.2 M ⊙ (800 pc). 2  Last, we model Spitzer 's limited angular resolution by removing stars in close proximity.When a source and its nearest neighbor (YSO or AGN) are within the adopted beam size threshold of 5 ′′ , we remove the lowermass source.We assign AGN a mass of 1.1 M ⊙ to avoid losing them to the mass cutoff.We do only one pass to remove sources, but this is sufficient to remove the vast majority of close neighbors.

Gas
We construct 2D projected column density maps with cloud distances of 200, 400, and 800 pc, which are chosen to model the majority of the clouds in the P20 sample.Figure 1 shows one of these maps with a spatial distribution plot of YSOs and AGN contaminants.
The Spitzer and Herschel fields of view focus on regions of high column density (clumps) within the clouds.To simulate this, we crop the gas maps to the bounds  set by a 10 21 cm −2 column density contour on a 120 ′′smoothed gas map constructed specifically for this purpose.This map is not used again in the further analysis.We smooth to keep small overdensities from ar-tificially enlarging the cropping area.3This greatly reduces the field of view compared to the full view, as shown in Figure 2, and makes our maps more similar to the Spitzer and Herschel data we compare with.Ad-ditionally, this significantly reduces the amount of lowdensity AGN contamination (see Figure 1 and Section 4 for more details).In order to simulate the angular resolution of Herschel, the gas maps are smoothed with a 36 ′′ Gaussian kernel.

Overview Statistics
To better compare with observations, We first define a few bulk cloud properties.We define the total cloud gas mass, M gas , as the combined mass of gas at column densities 10 21 cm −2 and above.Similarly, the dense gas mass, M dense , is the total gas at column densities 10 22 cm −2 and above.The dense gas mass fraction is then the ratio of dense to cloud gas mass.This metric gives an indication of the fraction of the cloud that is most likely to form clusters (Battisti & Heyer 2014;Heyer et al. 2016).We define the disk fraction as the ratio of the number of Class I and II YSOs to the total number of YSOs regardless of circumstellar material.Disk fraction can be used as a proxy for the population age (Haisch et al. 2001;Hernández et al. 2007b).A similar statistic, the Class II to Class I ratio, is generally believed to be a good relative evolution indicator for YSOs, especially for earlier evolution (G11, P20).
Figure 3 shows the evolution of the gas properties with time for the fiducial case.The cloud mass and dense gas mass increase steadily over time and peak at ∼ 5.4 Myr.The maximum mass reaches about half of the 20,000 M ⊙ of gas that makes up the entire simulated GMC.After this point, the cloud mass decreases rapidly to less than 1/10 of the initial GMC mass.The dense gas mass fraction exhibits a similar trend, peaking at M dense /M cloud ∼ 0.6 at the same time.
Figure 4 shows that the number of Class I and II sources evolve in a similar way to the gas.Star formation increases steadily for the first 3.43 Myr as indicated by the rising number of Class Is.After 3.43 Myr, star formation declines to 63% of the maximum in 2.45 Myr and then drops to only half this value in the next snapshot.The number of Class IIs evolves more gradually, peaking at 5.88 Myr, after which point it steadily decreases to about half its maximum value by the end of the simulation.formation event, which supports the use of a clusterderived model.It is shown with a vertical stretch and time axis shift to make the model more visible without adjusting the main model parameters.
We also plot a tweaked version of the model with minor adjustments to better fit our assumptions and outputs.Namely, we shift the time axis by 1.47 Myr to match our snapshots, increase the SFR from 100 to 435 (unitless metric which changes the vertical scale of the model), lengthen the rise and decay times for the Class Is from 0.5 to 1.7 Myr and 0.5 to 1.5 Myr, respectively, and shorten the lifetimes of Class Is and IIs to be closer to (but not exactly the same as) the half-lives for our adopted Class transitions (0.5 to 0.3 Myr and 2.0 to 1.5 Myr, respectively).With these parameters, the model reproduces the fiducial starforge data re-markably well.This suggests the starforge simulation provides a good representation of cluster formation.As we shall see below, the simulation appears to agree less well with star formation observed in full GMCs, which generally contain multiple distinct star-forming regions and have longer and more complex star formation histories.
Figure 5 shows the evolution of the Class II:I ratio and the disk fraction in this simulation.The disk fraction starts near 1.0 and then decreases nearly linearly to 0.21 in the final snapshot.This is more drawn-out than the traditional disk fraction versus stellar age plot (e.g.Mamajek 2009).The starforge calculation exhibits a broad range of Class II:I ratios, which span 1.3-19.0.For comparison, P20 recorded the Class II:I ratio and the cloud mass for 12 clouds at distances between 140 and 1400 pc.They found that the Class II:I ratio remained between ∼3.5-9.7 for each cloud observed, which is a much narrower range than we find in the starforge snapshots.However, the P20 values are uncorrected for AGN and edge-on disk contamination, which would likely change the Class II:I ratios, as will be seen below.
Using publicly available Herschel data (André et al. 2010, P20), we calculate the dense gas mass fraction of the clouds and clusters P20 and Gutermuth et al. (2009) observed.We adopt the publicly available YSO lists from SESNA, correct for AGN and edge-on disk contamination, and crop for coverage consistency and to the N(H 2 ) = 10 21 cm −2 limit.In the case of the Gutermuth et al. (2009) data shown, we adopt all "cluster cores" that overlap with clouds from the P20 sample, and crop to square areas that are twice the diameter implied by the R circ radii listed in that paper, once converted to the most recent heliocentric distances reported in P20.Some of the selected areas of adjacent cluster cores overlap significantly.The assumed and computed data for these plots are listed in Tables 3 & 4.
Figure 6a shows that starforge and the clouds in P20 occupy different regions of the Dense Gas Mass Fraction -Class II:I ratio parameter space.The trajectory agrees better with the clusters from Gutermuth et al. (2009), except for the earliest snapshots.We could correct for this by assuming some amount of ambient star formation occurs in the cloud before the main cluster forms, which would increase the Class II:I ratio, more noticeably in the early and late snapshots that have few Class I and II sources.This supports the implication that starforge more closely models the formation of a large cluster rather than star formation in a GMC.Inspection of Tables 3 & 4 indicates that the total gas mass and dense gas mass in the simulation are also more consistent with the ranges reported for the Gutermuth et al. (2009) clusters.We next apply a correction for AGN to the synthetic observations by removing 4.5 sources per square degree for both Class Is and IIs.We find that the synthetic observation trajectories exhibit strong agreement with each other and the fiducial (Figure 6b).This is expected, since we add that same density of AGN contaminants at the beginning of the synthetic analysis.

Evolution of the Star-Gas Fraction
The calculation of the S-G correlation in this work emulates the treatment from P20, allowing us to better compare the outcomes of the two.We calculate the n th nearest neighbor distance (NND) for each Class I YSO, for n = 11 using scipy.spatial.KDTree.KDTree uses the algorithm described by Maneewongvatana & Mount (1999) to create a binary tree of 3-dimensional nodes (the positions of the sources).This allows for the quick lookup and classification of nearest-neighbors.We use n = 11 because it is a good compromise of spatial resolution (typically 0.1-2 pc smoothing-scale in nearby clouds) and low relative uncertainty (33%, Casertano & Hut 1985).This choice is consistent with Casertano & Hut (1985), G11, and P20.Using a circular mask with a radius of the NND we calculate the area A of each circular mask, the mean column density in each circle Σ C , and the ratio C of covered area to total area within each circle.C covers edge effects and is thus almost always unity.From this, we calculate Σ gas , the gas mass surface density.Σ YSO is the measurement of the surface density of YSOs, calculated as (Casertano & Hut 1985) where M YSO is the adopted mean mass per YSO and n, A n and C are defined above.Except for our fiducial analysis, where we try to avoid as much observational bias as possible, we fix M YSO at 0.5 M ⊙ to keep the analysis consistent with P20.
Figure 7 shows the median Σ YSO /Σ 2 gas versus time, which captures the vertical offset and spread around the power-law fit (G11).While anything more than a general positive trend with time showing increasing stellar density as a function of gas density is not immediately clear, the Class I and II values are close to each other for the first ∼ 6 Myr.After this point, the populations no longer appear correlated.This points to a large-scale decoupling of the YSOs from their surrounding gas at around 6-6.5 Myr.This is supported by visual examination of the snapshots.Figure 8 shows a snapshot before decoupling occurs (3.42 Myr) and a snapshot after decoupling occurs (7.82 Myr).It is clear that nearly all YSOs reside near or within dense gas before the decoupling, but afterwards the two groups are significantly less correlated.
The dense gas mass fraction peaks at ∼5.38 Myr (Figure 3).This is when feedback begins to disperse the cloud (see Figure 2), and there is a ∼1 Myr lag before the effects are seen in the other statistics.For example, Figure 5 shows that the number of Class Is drastically declines and the number of Class IIs peaks at ∼6.36 Myr.This causes the Class II:I ratio to rise significantly.And, as mentioned above, this is also the time when Class Is and IIs in Figure 7 appear to decouple.

Star-Gas Correlation versus Time
Figure 9 shows the slopes and uncertainties of the S-G correlations for the fiducial analysis along with the three sets of synthetic observations as a function of time.Most of the slopes lie relatively close to 2.0, however the wellcorrelated slopes either lie above or below 2.0, usually localized around 2.4-2.5 or 1.7-1.8.Over half of the fiducial snapshots visually appear to have a tight YSO and gas surface density correlation, with an uncertainty in the slope of ≤ 0.2.This provides significant evidence that the power-law relationship for the S-G correlation is a real effect that is a result of underlying physics and not a result of observational bias (see Figure 13 in the Appendix for the fiducial S-G correlation plots).
However, many of the snapshots are not wellcorrelated, appearing as a clump of points that lie on the for AGN and edge-on disk contamination, and cropped for coverage consistency and N(H2) > 10 21 cm −2 , so they differ to varying degrees from the raw values reported in those works.Black points and line represent the time evolution trajectory for the fiducial analysis of this work, starting from the bottom left.Right: fiducial trajectory overlaid with trajectories from the synthetic analyses at different distances, corrected for AGN contamination.Note that points at high Class II:I ratio are highly uncertain (e.g. Figure 5c). .Median value of ΣYSO/Σ 2 gas versus time for both Class I and Class II sources.In addition to showing the increasing stellar density as a function of gas density, these values are closely correlated until ∼6 Myr (dashed vertical line).At this time, feedback clearly begins to disrupt the gas (see Figure 3) thereby inducing decoupling of the gas structure and the YSO distribution.
expected line, but do not span a significant range of surface densities.This is especially true for snapshots with fewer than ∼100 Class I sources, since this often leads to poorly-constrained slopes with error bars as large as 0.6.This difficulty hinders comparison with previous observations, as many of the observed clouds in G11 and P20 have many more sources and more completely populate the S-G space and thus the S-G correlation.However, the addition of synthetic observation effects, especially adding AGN or removing close neighbors, can artificially compensate for this by filling out the low-density region and depleting the high-density region of the plot, respectively.This is discussed in Section 3.4 below.In addition, the S-G correlations for each snapshot of the simulation in the fiducial and 200, 400, and 800 pc synthetic analyses can be found in the Appendix.
Figure 9 illustrates the evolution of the S-G slope as a function of time in the simulation.While the shape, slope, or scatter of the S-G correlation do not change monotonically with time, there are several features that are roughly independent of distance and the presence of the synthetic considerations.This implies that the synthetic observation effects don't obscure the underlying physics, except for in snapshots where low-number (of YSOs) statistics are significant (e.g., the poorly correlated snapshot in Figure 12).Figure 9 shows the S-G slope declines until around ∼ 4 Myr (most noticeable at closer distances), at which point the number of Class I sources peaks.From ∼4-6 Myr the S-G slope increases as the number of Class II sources continues to rise; the peak in the S-G slope at 6 Myr coincides with the peak in the number of Class II sources.After 6 Myr the S-G slope declines sharply as feedback begins to disperse the cloud in earnest.Many of the snapshots around this time also have poor S-G correlations.Even though much of the cloud gas is dispersed from the central region, the YSOs' dynamics take longer than the gas to respond to the changing gravitational potential.However, star formation still occurs in the remaining pockets of dense gas, maintaining some degree of S-G correlation in the later snapshots (Figure 1).While it is clear that some star-gas statistics evolve with time, the slope of the S-G correlation appears to be relatively constant across the history of the cloud.This is consistent with the observations of G11 and P20 who found little variation in the slopes across a wide collection of GMCs with very different ages.The spread in the S-G slopes and fit uncertainties are significantly larger for the simulations than for the observations of G11 and P20, however.This comparison may not be fully equal, as there is a selection effect on which clouds are actually observed and included for analysis.For this reason the observed clouds may span a narrower range in the cloud lifetime: young clouds with little star formation may not be identified as distinct and/or interesting star-forming regions and thus will be excluded, while older clouds that are in the process of dispersing are excluded since they have little remaining dense material.As a result, the especially-poorly-populated early and late snapshots are not well represented in real data, as it is difficult to find and observe very young and very old clouds.More observational work will need to be done to more effectively compare with these snapshots.

Demonstration of Synthetic Effects on the Star-Gas Correlation
In this section, we explore how each synthetic effect impacts the apparent S-G correlation.Figures 11 and 12 compare the fiducial S-G correlation with those obtained for five different synthetic effects.
The first effect we add to the synthetic observation is the adoption of a uniform YSO mass. Figure 10 shows the mean and spread of YSO masses in the simulation, and as can clearly be seen, adopting a fixed average mass does not well-represent the true average mass, which varies by a factor of ∼ 10 over time.However, since in- In many of the earlier snapshots, undersampling causes large uncertainties and produces slopes that are discrepant with observations.We limit the y axis to better compare differences between the runs.This obscures the (unreasonable) points of some of the earlier snapshots.See Figures 13-16 for full slope and uncertainty values.The leftmost point in the bottom right panel is missing since there is no slope for that individual snapshot (see Figure 16).
dividual real YSO masses cannot be directly measured, observational analyses such as P20 must adopt some approximation.Figures 11b and 12b illustrate the S-G correlation assuming a uniform YSO mass of 0.5 M ⊙ .
Comparing panels (a) and (b) indicates that using the true mass of the sources has little effect on the S-G correlation.While the points move slightly vertically, the slopes change by less than 0.1.Consequently, source mass appears to have a relatively minor impact on the S-G correlation.Considering the minor effects, the uniform mass is used in the rest of the demonstration.The impact of the removal of close neighbors is more significant.When multiple close sources appear as a single source, the effect is to remove many of the highest density points in the S-G relationship shown in Figure 9c.This, in turn, has a flattening effect on the power-law slope, for example, bringing the slopes of most snapshots (all snapshots between 2 and 6.5 Myr) within 2σ of 2.0 (see Figure 9).The earlier and later snapshots tend to be (often significantly) less well-sampled, which likely explains their inconsistent slopes (as shown in Figure 12).The impact of this effect increases with distance, as the 5 ′′ minimum separation imposed on the YSO lists translates to larger physical separations.This can be observed qualitatively in Figures 13 -16 in the Appendix.
Figure 11d illustrates the effect of detection limits on the S-G correlation.We find that implementing a mass sensitivity limit significantly decreases the number of sources at all densities, which increases the fit uncertainties across all snapshots.However, this does not significantly change the slope in well-correlated snapshots.In contrast, Figure 12 shows that the addition or removal of a single point can significantly change the slope in a snapshot with fewer YSOs.This effect is also more extreme at larger distances.Specifically, at 200 and 400 pc, the mass limit is the Hydrogen-burning limit: 0.08 M ⊙ , while the limit at 800 pc is 0.2 M ⊙ , significantly reducing the number of sources with which to calculate the S-G correlation.Compare Figure 16 with Figures 14 and 15 in the Appendix for a visualization of how the number of sources decreases with increasing distance.
Next we investigate the impact of AGN contamination on the S-G correlation.The addition of AGN has a significant impact on the S-G correlation as shown in Figures 11e and 12e.Since the AGN are uniformly randomly distributed throughout the field of view, they add a relatively constant (Σ YSO = ∼0.3,∼0.1, ∼0.03 M ⊙ /pc 2 at 200, 400, and 800 pc, respectively) "foot" of points to the bottom of the S-G correlation.This disproportionately affects low Σ YSO regions and artificially flattens the power law.The flattening increases with distance, so much so that the slope of the S-G correlation for every snapshot at 200 pc and 400 pc would lie below 2 at all times.However, we follow observational convention and implement a column density cutoff when fitting the slope as described below.
For nearby clouds, the number of YSOs observed at relatively low column densities is small, which causes observations of those regions to be dominated by AGN contaminants.To deal with the similar issue of our syn-thetic AGN, we adopt the same approach as G11 and remove YSOs in our catalog in regions with log(Σ gas ) < 1.3 M ⊙ pc −2 and refit the remainder.This is demonstrated in Figures 11e,f and 12e,f.Applying this treatment to the synthetic observations with AGN confirms that such a cut is justified to minimize the bias of the fit caused by AGN contamination.After applying this cut, most slopes steepen and approach the expected value of ∼2.0 (Figure 9).The presence of many well-correlated S-G relationships for the fiducial starforge run implies that the S-G correlation is a physical phenomenon and not solely the result of observational biases.However, the addition of synthetic observation considerations does artificially lower the slope of the S-G correlation for many of the snapshots, generally increasing agreement with observations.This begs the question of whether the very consistent value of 2 determined by P20 is partially caused by observational effects.In this case, the S-G correlation slope is not as invariant and universal as it appears in P20.However, in the latest snapshots, which have better agreement in Class II:I ratio versus dense gas mass fraction with P20 clouds (Figure 6), the S-G correlation slopes are much lower than observed.
Nonetheless, it is striking that the broad range of evolutionary stages spanned by one cloud, modeled with all key physical effects, produce a relatively uniform powerlaw slope.Once star formation is underway, stellar feedback helps to regulate the relationship between dense gas and YSOs.Clouds with particularly high Class II to I ratios are likely dominated by stellar feedback and in the process of cloud dispersal.Follow-up observations that minimize observational effects are required to fully constrain if and the extent to which these biases conspire to produce an S-G power-law slope of ∼ 2.

Comparison to Previous Work
Chevance et al. ( 2022) find that GMCs in nine nearby disc galaxies usually disperse within ∼ 3 Myr after unembedded high-mass stars emerge.While not directly measured in this work, we believe the first high-mass stars likely emerge shortly before feedback begins dispersing the cloud in earnest.We estimate dispersal to become qualitatively significant sometime between ∼ 5.4 − 6.4 Myr, as described in Section 3. And, considering the GMC is nearly completely dispersed by 8.8 Myr, the simulations are consistent with the observed ≲ 3 Myr time frame, as well as the ∼ 10 Myr total cloud lifetime they estimate.There is a slight vertical shift in the points on the plot, but it does not significantly change the slope or the uncertainty of the S-G correlation.We adopt uniform mass for the rest of the panels.c) The removal of close neighbors that would have been indistinguishable by Spitzer.This predominantly removes high-density points, lowering the slope.d) The removal of low-mass sources that would have been undetected by Spitzer.This removes points without visible bias towards density, increasing the uncertainty in the slope.e) The addition of AGN.This predominantly adds low-density sources, lowering the slope.f) All previous synthetic effects at once.The slope is much closer to 2. Black dashed lines represent a density cutoff imposed to account for the presence of AGN.The slopes and best-fit lines for e) and f) are only based on points to the right of the black line.
As mentioned in the Introduction, P20 adapted HD simulations by Qian et al. (2015) to create synthetic observations of their 12 observed clouds.These synthetic observations included 2D projection and neighbor removal.The HD simulations produced slopes between 2.3 − 2.7, higher than the observed 1.8 − 2.3.The simulations are also limited in density dynamic range compared to some of the clouds they model (see Figure 6 in P20).However, the simulated slopes are similar to the values of 2.0 − 3.0 we observe in the fiducial run (before cloud dispersal and excluding the first two snapshots, see 2).Caution is required when comparing with these simulations since P20 modeled 12 different clouds at a single time, while this work models one cloud at many different times.Regardless, the main improvement of starforge over the simulations in Qian et al. (2015) is more realistic physics, especially magnetic fields and kinematic feedback.While magnetic fields do not play a very significant role in setting the slope of the S-G correlation, kinematic feedback allows starforge to evolve the GMC without driven turbulence (which was necessary for the simulations in Qian et al. 2015).While the STARFORGE simulation starts with an initial turbulent setup, the turbulence, evolution, and dispersal of the cloud are regulated entirely by stellar feedback.

Model and Analysis Caveats
While starforge faithfully reproduces the S-G correlation in many snapshots, there are some areas for improvement.For example, Figure 6, which displays the Class II to I ratio versus dense gas mass fraction, that the simulation exhibits poor agreement with the P20 clouds.In contrast, the simulation agrees better with the cluster data from Gutermuth et al. (2009).This implies that the starforge simulation analyzed here produces something that is closer to a smaller denser structure (i.e., a cluster) than the stellar complexes formed in the GMCs of the P20 sample, which may be characterized by a longer and richer star formation history.
The drastic evolution in the number of Class Is with time highlights that the simulation SFR is not constant, in contrast with the assumption of constant SFR made by G11, P20, and others when using class ratios (to infer age, for example).Figure 5 shows that this leads to Class II:I ratios that vary much more than in the P20 observations.Megeath et al. (2022) argued that a variable SFR similar to that produced by the starforge simulation is necessary to explain the ensemble of Class II:I ratios and disk fractions in nearby clusters.The agreement between this model and the starforge data (Figure 4) provides more evidence that starforge produces something more similar to a monolithic cluster than a full GMC (i.e., with several smaller, distinct clusters).
However, we caution that here we only analyze one simulation that aims to model the typical conditions of a Milky Way cloud.Future work is needed to explore the broad range of conditions modeled in the starforge simulation suite, which includes clouds with varying initial magnetic field, turbulence, interstellar radiation field, surface density, cloud size, and cloud initialization (Guszejnov et al. 2022).In particular, the initial cloud setup, that of a uniform density sphere, is a significant oversimplification of the complexity of forming and accreting molecular clouds.Overall, agreement with both data sets would likely be improved by using more realistic initial conditions.For example, a slower star formation start could increase the Class II:I ratio in the early snapshots, improving agreement.Simulations that begin with more realistic cloud initialization, such as a driven-turbulence sphere (Lane et al. 2022) or models that follow cloud formation from galaxy scales (Hu et al. 2023;Ganguly et al. 2023, Hopkins et al. 2023 in prep.), are likely necessary to advance agreement between the starforge framework and observations.
One recent interesting aspect of starforge comes from Grudić et al. (2023), who ran 100 2000 M ⊙ STAR-FORGE simulations and found a sharp mass cutoff on the IMF at 28 M ⊙ .This is in contrast to a simulation with similar parameters but 10 times the mass, which generated a 44 M ⊙ star, and a simulation with 10 times the gas surface density which generated a 107 M ⊙ star.They suggest that the STARFORGE IMF has a highmass cutoff that depends on the environment.This cutoff is generally different from the canonical 100−150M ⊙ cutoff, which they conclude implies that the IMF cannot be reproduced in small clouds simply by randomly sampling from the full IMF.
Here we also outline some inconsistencies in our data processing.The bounds on the cropped field of view (see Figure 1) are set by the furthest extent of the N(H 2 ) = 10 21 cm −2 contour.This occasionally causes larger-than-intended fields of view when an area of gas denser than N(H 2 ) = 10 21 cm −2 is present away from the central cluster.The only major impact this has on the analysis is to increase the number of AGN when calculating the S-G correlation.However, this impact is largely mitigated by the low density cut discussed previously.We also neglect a number of steps that would be needed to complete a fully "apples-to-apples" comparison with the observations.For example, we do not use radiative post-processing to construct the YSO SEDs (e.g., Offner et al. 2012).Nor do we construct synthetic dust continuum maps in order to compute the column density (e.g., Juvela 2019).These steps would allow us to apply the observational biases, such as the detection limits, more directly.However, we expect any impact on the S-G slope to be minimal, since the YSO positions and relative amounts of dense gas would be unchanged.

CONCLUSIONS
In this study, we examine a 20,000 M ⊙ star-forming cloud in the starforge simulation suite in order to investigate the presence and evolution of the S-G correlation.To effectively do so, we create synthetic observations to compare with previous observational work, specifically P20 and G11.These synthetic observations include 2D projection of gas and star particle distributions at multiple distances, an age-to-Class conversion for the simulated stars using an exponential decay model, AGN contamination, a low stellar mass cutoff, the removal of close (unresolved) neighbors, gas map smoothing to mimic limited angular resolution, and a field-of-view crop at a gas column density of N(H 2 ) = 10 21 cm −2 .
Since most of these effects depend on distance, we place each cloud at 200, 400 and 800 pc to mimic the distances of star-forming regions observed by G11 and P20.This changes the angular size of the cloud, the number of AGN, the mass sensitivity limit, and the neighbor threshold.From these synthetic observations, we examine the dense gas fraction, YSO distribution and frequency, and the S-G correlation for the fiducial analysis and for the synthetic analyses at each distance.
We find that the starforge simulation successfully reproduces the S-G correlation in many snapshots and exhibits a typical S-G slope within 1σ of the observed slope of 2. The presence of the S-G correlation both with and without accounting for observational effects implies that this is a real relationship that is a product of the underlying physical processes.However, observational biases, such as AGN contamination, appear to strengthen the S-G correlation, reduce time variation and promote a slope closer to 2.
We find that the Class II:I ratios and dense gas fraction characteristic of the starforge simulation exhibit better agreement with those of the clusters in the Gutermuth et al. (2009) sample than the stellar complexes forming in the clouds in P20.No regions in either observational study match the low Class II:I ratios found at early times (< 3 Myr) in the simulation.This implies that the P20 and Gutermuth et al. (2009) clouds/clusters form stars at a low rate for a few million years.Thus, bias in cloud selection, which favors actively star-forming clouds with significant amounts of dense gas, possibly also contributes to the apparent universality of the S-G correlation.
The present study only considers the S-G correlation under one set of typical simulated cloud conditions.Future work is needed to examine the impact of cloud properties and more realistic initial conditions on the S-G correlation.

A. TABLES
In Tables 1 and 2, we present various statistics extracted from the analysis of the fiducial snapshot.We present statistics of the clouds and cluster cores from P20 and G09 used in our analysis in Tables 3 and 4. We have updated the data from P20 and G09 to the latest datasets from SESNA and Herschel, corrected for AGN and edge-on disk contamination, and adopted distances from P20.

B. S-G CORRELATION PLOTS
We present here in Figures 13,14,15,and 16 the full collection of S-G correlation plots for the fiducial and synthetic analyses.The fiducial case represents analysis done with minimal adjustments, while the others contain all synthetic effects described in Section 2.2.

Figure 1 .
Figure 1.Projected N(H2) column density map of a 200 pc-distance cloud with N(H2) = 10 21 cm −2 contour over-plotted in green.Colored circles indicate the locations of YSOs and AGN.a) Full field of view of the simulation at ∼5.4 Myr.b) Zoomed (∼20-pc) field of view cropped to the furthest extent of the green contour at ∼5.4 Myr.AGN contaminants dominate the source counts in the low-column density regions.

Figure 2 .
Figure 2. Projected N(H2) column density map with N(H2) = 10 21 cm −2 contour over-plotted in green and N(H2) = 10 22 cm −2 contour in magenta at (a) ∼2.4 Myr.(b) ∼5.4 Myr.(c) ∼8.3 Myr.The green contour outlines the likely Spitzer field of view for an equivalent cloud.Note that the high-density (magenta) region coalesces as star formation increases and eventually breaks apart due to stellar feedback, which is in the process of dispersing the cloud in (c).

Figure 4 .
Figure 4. Evolution of each Class of starforge-derived YSO counts in this work (black points) overlaid with analytical models adapted from Megeath et al. (2022): (a) Number of Class I sources versus time, (b) Number of Class II sources versus time, (c) Number of Class III sources versus time.The orange lines are shifted and rescaled versions of the Megeath et al. (2022) models using their parameter selections, while the blue lines adopt parameter value adjustments to achieve strong agreement with the STARFORGE data.

Figure 5 .
Figure 5. (a) Class II:I ratio versus time increases steadily until about 6 Myr, at which point it jumps up and does not follow a consistent trend.(d) Disk fraction, i.e., ratio of Class I and Class II sources to total number of sources, decreases steadily, but slower than comparable observations based on mean stellar ages in real clouds.Error bars are calculated through standard error propagation.

Figure 6 .
Figure 6.Dense Gas Mass Fraction versus Class II:I ratio.Left: Blue triangles are values from nearby molecular clouds in P20.Orange triangles represent clusters in those clouds fromGutermuth et al. (2009).All observed data have been corrected for AGN and edge-on disk contamination, and cropped for coverage consistency and N(H2) > 10 21 cm −2 , so they differ to varying degrees from the raw values reported in those works.Black points and line represent the time evolution trajectory for the fiducial analysis of this work, starting from the bottom left.Right: fiducial trajectory overlaid with trajectories from the synthetic analyses at different distances, corrected for AGN contamination.Note that points at high Class II:I ratio are highly uncertain (e.g.Figure5c).
Figure 7. Median value of ΣYSO/Σ 2gas versus time for both Class I and Class II sources.In addition to showing the increasing stellar density as a function of gas density, these values are closely correlated until ∼6 Myr (dashed vertical line).At this time, feedback clearly begins to disrupt the gas (see Figure3) thereby inducing decoupling of the gas structure and the YSO distribution.

Figure 8 .
Figure 8. Projected N(H2) column density map with N(H2) = 10 21 cm −2 contour over-plotted in green and N(H2) = 10 22 cm −2 contour in cyan at (a) ∼3.4 Myr.(b) ∼7.8 Myr.Note in (a) that most YSOs, especially Class I sources, remain close to or within the dense gas cyan contour.In (b) however, many of the YSOs are no longer correlated with the locations of the denser gas at either contour level, indicating YSO-gas decoupling.Existing YSOs (mainly Class IIs and IIIs) remain relatively stationary for the first few million years as the gas dissipates, being bound together by gravity.However, new Class Is continue to form in the denser gas (almost all Class Is are within the cyan contours).

Figure 9 .
Figure9.Slope and uncertainty of the S-G correlation for each snapshot.The shaded region corresponds to the range of values observed by P20.In many of the earlier snapshots, undersampling causes large uncertainties and produces slopes that are discrepant with observations.We limit the y axis to better compare differences between the runs.This obscures the (unreasonable) points of some of the earlier snapshots.See Figures13-16for full slope and uncertainty values.The leftmost point in the bottom right panel is missing since there is no slope for that individual snapshot (see Figure16).

Figure 10 .
Figure10.Average combined Class I and II mass for each snapshot.Error bars represent 95th percentile.It is clear that an assumed mass of 0.5 M⊙ does not accurately represent YSO masses at all times.However, this has little effect on the calculation of the S-G correlation (see Section 3).
Implications for the S-G Correlation

Figure 11 .
Figure11.Comparison of different synthetic observational effects on a well-correlated snapshot.a) The "fiducial" analysis with no extra considerations.Each panel b) through e) demonstrates 1 synthetic effect each.b) In the calculation of the S-G correlation, uniform 0.5 M⊙ mass for each source is used.There is a slight vertical shift in the points on the plot, but it does not significantly change the slope or the uncertainty of the S-G correlation.We adopt uniform mass for the rest of the panels.c) The removal of close neighbors that would have been indistinguishable by Spitzer.This predominantly removes high-density points, lowering the slope.d) The removal of low-mass sources that would have been undetected by Spitzer.This removes points without visible bias towards density, increasing the uncertainty in the slope.e) The addition of AGN.This predominantly adds low-density sources, lowering the slope.f) All previous synthetic effects at once.The slope is much closer to 2. Black dashed lines represent a density cutoff imposed to account for the presence of AGN.The slopes and best-fit lines for e) and f) are only based on points to the right of the black line.

Figure 12 .
Figure 12.Comparison of different synthetic observational effects on a poorly-correlated snapshot.Features of figure the same as in Figure 11.Note that this snapshot is particularly sensitive to the removal of a single high-density point.

a
Stellar completeness corrected by 0.163 following P20.

a
Stellar completeness corrected by 0.163 following P20.

)Figure 14 .Figure 15 .Figure 16 .
Figure13.S-G Correlation Plots for each snapshot in the fiducial analysis.Note that the first snapshot is extremely undersampled.

Table 1 .
Table of various fiducial snapshot statistics.