Testing the Relationship between Bursty Star Formation and Size Fluctuations of Local Dwarf Galaxies

Flat cosmological models with a mixture of dark energy, cold dark matter, and baryonic matter are consistent with the observations of structure formation on large scales ( ≫ 1 Mpc; Bahcall et al. 1999; Frenk & White 2012; Vogelsberger et al. 2014; Schaye et al. 2015; Planck Collaboration et al. 2020). In this cold dark matter model (ΛCDM), dark matter is viewed as particles with weak self-interactions, as well as negligible interactions with baryonic matter. However, on galactic scales (<1 Mpc) there are several tensions between observations and predictions of ΛCDM (for a comprehensive review see Bullock & Boylan-Kolchin 2017).

Dwarf galaxies, with stellar masses ≲ 10⁹ M_⊙, are known to challenge CDM predictions on small scales (Brooks 2019). It has long been recognized that the observed dark matter halos of dwarf galaxies exhibit a different core density profile (Flores & Primack 1994; Moore 1994) than the predicted dark-matter-only Navarro–Frenk–White (NFW; Navarro et al. 1997) profile.

The NFW density profile rises steeply at small radii (Navarro et al. 2010). However, observations of rotation curves of dwarf galaxies reveal rather a constant-density profile in apparent contrast with the NFW profile (McGaugh et al. 2001; Marchesini et al. 2002; Simon et al. 2005; Kuzio de Naray et al. 2008; Oh et al. 2015; Relatores et al. 2019), referred to as the "cusp–core" problem. Two solutions have been proposed in the literature.

The first solution is a modification to the ΛCDM in which dark matter particles are posited to be strongly self-interacting (Spergel & Steinhardt 2000; Kaplinghat et al. 2016). In the self-interacting dark matter model, dark matter collisions allow energy exchange between particles and thermalize the inner halo over the cosmological timescale, leading to a cored central density (Vogelsberger et al. 2012; Peter et al. 2013; Elbert et al. 2015; Fry et al. 2015; Sameie et al. 2020).

The second solution to the cusp–core problem is to implement baryonic feedback (in forms of supernova explosions, radiation pressure, photoelectric heating, photoionization, and stellar winds; Navarro et al. 1996; Governato et al. 2010; Sales et al. 2014; Hopkins et al. 2020). In these models, outflows that are driven by baryonic feedback over the course of star formation (SF) are capable of inducing fluctuations in the central gravitational potential of the galaxy, transferring energy into dark matter particles, and ultimately flattening the dark matter core density profile in dwarf galaxies (Governato et al. 2010; Sales et al. 2010; Pontzen & Governato 2012; Chan et al. 2015; Read et al. 2016; Tollet et al. 2016). This alternative solution requires high-resolution hydrodynamical simulations to model the inhomogeneous interstellar medium and couple the SF to high-density gas regions (Hopkins et al. 2014, 2018). Implementation of baryonic feedback into CDM-only simulations has been shown to further explain several other characteristics of dwarf galaxies such as their low stellar-to-halo mass ratio (Hirschmann et al. 2012; Behroozi et al. 2013; Moster et al. 2013).

Other than forming a cored dark matter density profile, baryonic feedback is predicted to have additional impacts on the dynamical and morphological properties of dwarf galaxies (El-Badry et al. 2016, 2017; Read et al. 2019), many of which are shown to be in agreement with recent observations. For instance, outflows caused by supernova explosions result in bursty SF rates (SFRs) in dwarf galaxies (McQuinn et al. 2010; Sparre et al. 2017), which is consistent with the different SFRs inferred from SFR indicators that probe different timescales in local and high-redshift dwarf galaxies (Weisz et al. 2012; Domínguez et al. 2015; Guo et al. 2016; Mehta et al. 2017; Caplar & Tacchella 2019; Emami et al. 2019; Faisst et al. 2019; Flores Velázquez et al. 2021). Outflows can also explain the lack of coherent disk formation in most irregular dwarf galaxies (Wheeler et al. 2017), as well as an observed specific SFR–gas velocity dispersion (sSFR–σ) correlation in this mass range (Stilp et al. 2013; Cicone et al. 2016; Hirtenstein et al. 2019; Pelliccia et al. 2020). Furthermore, outflowing winds that are powered by episodic bursts of SF drive radial migration of star and gas particles over the course of SF (De Young & Heckman 1994; Graus et al. 2019; Mercado et al. 2021), which is consistent with metallicity and inverted age gradients in low-mass galaxies (Aparicio & Tikhonov 2000; Stinson et al. 2009; Hidalgo et al. 2013; Vargas et al. 2014; McQuinn et al. 2017; Jafariyazani et al. 2019).

The predicted radial migration has another consequence: the outflowing winds can also impact the size of the galaxy. In particular, during bursts of SF gas and stars move into the center of the galaxy until the first supernova explodes. Subsequent supernova explosions will yield a shallow gravitational potential, gas and stellar migration toward larger radii resulting in more extended galaxy morphologies. Over time, the expelled gas will cool and collapse back onto the galaxy center, triggering the next burst of SF. This causes size fluctuations during bursts of SF across the galaxy (El-Badry et al. 2016). Yet this prediction has not been observationally confirmed for the real dwarf galaxies.

In this work, we aim to test this prediction for a sample of local dwarf galaxies with a stellar mass range of 10⁷–10^9.6 M_⊙, the same regime suggested by El-Badry et al. (2016). We will study the relationship between bursty SF and size fluctuations in both observed and simulated galaxies and provide physical intuitions for our findings.

An overview of this paper is as follows. We describe the observed sample and our measurements of the galaxies' stellar mass, dust-corrected luminosities, and R-band size in Section 2. In Section 3 we discuss our observed findings. In Section 4 we analyze a suite of simulated dwarf galaxies, analogous to our observed sample, from FIRE-2 simulations and provide physical interpretations about the effect of bursty SF on the galaxy's size fluctuation. We discuss previous works and the implications for the era of the James Webb Space Telescope in Section 5. Lastly, we summarize our results in Section 6.

Throughout this work we assume that the default ΛCDM cosmology has parameters H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_m =0.315, and Ω_Λ = 0.685 (Planck Collaboration et al. 2020).

In this work, we use a sample of galaxies from the Local Volume Legacy (LVL) survey (Kennicutt et al. 2008). LVL is an unbiased and statistically complete sample of 258 nearby field star-forming galaxies. Galaxies in the LVL sample are selected based on a set of criteria including galaxies that are within 11 Mpc (D ≤ 11), lie outside the Galactic plane (∣b∣ > 20°), are brighter than B = 15 mag, and span an RC3 type galaxy range T ≥ 0 (i.e., galaxies of spiral and irregular morphologies later than S0a). Due to its volume-limited nature, the LVL sample is dominated by dwarf galaxies and is complete above 10^6.5 M_⊙ stellar mass. The LVL public data set consists of GALEX ultraviolet (Lee et al. 2011), narrowband Hα imaging (Kennicutt et al. 2008), Spitzer IRAC and MIPS imaging (Dale et al. 2009), optical MMT spectroscopy (Berg et al. 2012), and optical ground-based imaging (Cook et al. 2014).

About 70% of the sample has deeper R-band imaging (∼0.5 mag deeper than the Sloan Digital Sky Survey (SDSS); Cook et al. 2014) obtained from 1–2 m telescopes located at Cerro Tololo Inter-American Observatory, Kitt Peak National Observatory, Vatican Advanced Technology Telescope, and the Wyoming Infra-Red Observatory. The median pixel scale of the instruments is 05 pixel⁻¹, with a standard deviation of 01 pixel⁻¹. For the other 30% of the sample, we use SDSS r-band imaging taken from SDSS DR7 (Abazajian et al. 2009). We use the data product produced in Cook et al. (2014). All the images are sky subtracted and corrected for the contaminant sources, e.g., background galaxies and foreground stars.

The Hα flux measurement is described extensively in Kennicutt et al. (2008). Different aperture configuration methods were applied on galaxies in the sample to ensure that the Hα flux is integrated over the entire galaxy, not a subset.

The stellar masses and dust extinctions are obtained from Weisz et al. (2012) and Johnson et al. (2013). Briefly, the masses are derived by fitting the observed UV, optical, and IR luminosities of each galaxy with the Bruzual & Charlot (2003) suite of stellar population synthesis models and a varying set of parameters (age, stellar metallicity, exponentially declining SF histories (SFHs), and dust). The Hα and UV dust extinctions are derived from the "energy balance" extinction correction method (Kennicutt et al. 2009; Hao et al. 2011), by which the intrinsic luminosity is retrieved from the combination of unobscured (e.g., Hα and UV) and obscured (far-IR continuum) signatures of SF. For the 30% of the sample where far-IR measurements are not available, a correlation-based attenuation correction is used (Lee et al. 2009). For that, the Hα dust attenuation is determined via a scaling relation between A_Hα and B-band magnitude (M_B ; Lee et al. 2009). As shown in Lee et al. (2009), the scatter between A_Hα and M_B decreases with decreasing luminosity, such that for the relatively transparent dwarf galaxies that dominate our sample (−18< M_B < −14.7) the scatter is roughly 10%–20%. The Hα extinction corrections used here yield similar results to those used by Lee et al. (2009). The dust corrections derived by the application of energy balance and the empirical A_Hα –M_B correlation are much smaller than the observed scatter in the Hα luminosity and do not have a significant effect on our interpretations about the bursty SF of the sample. The UV dust attenuation (A_UV) is also estimated as the Hα attenuation (A_Hα ) scaled by 1.8.

We determine the sizes of the galaxies in the R-band images. This is because the galactic potential fluctuations caused by the outflowing winds impact the distribution of all stars of any age. Additionally, through the R-band imaging, light from both young and old stellar populations can be assessed via the combination of strong nebular lines (i.e., Hα, which traces stars that are formed a few Myr ago) and the stellar optical continuum (at a wavelength range of 5800–7300 Å, which probes stars that are formed over hundreds of Myr). We run Source Extractor (Bertin & Arnouts 1996) to determine the half-light radii in the R band. Since the images are sky subtracted, we set the background value to zero in the code. For each galaxy, the zero-point, pixel scale, and gain parameters are set by values stored in the header files. These parameters depend on the characteristics of the telescope with which the object was observed. We set the FWHM parameter to the seeing of the observation. We choose a flexible elliptical aperture to measure the flux of the detected objects inside that aperture as described in Kron (1980). Whenever the image is highly resolved, which is mostly the case for objects that are closer than 3 Mpc, the Source Extractor algorithm segments the galaxy into multiple sources. In these cases we smooth the image by applying a Gaussian kernel with a minimum FWHM necessary for Source Extractor to recognize the galaxy as a single source. It may be of concern that the smoothing procedure washes out the galaxies' substructures, such as the spiral arms; therefore, the measured half-light radii are smaller than when the smoothing is not applied. However, the convolution that we perform is with a much (factor of ∼5) smaller profile than the ultimate sizes measured in these galaxies. This means that the measured sizes will not be significantly increased by performing this convolution. Source Extractor yields the half-light radius from a circular photometric aperture within which 50% of the total flux of the galaxy is enclosed. The choice of circular aperture is an acceptable approximation since the galaxies are mostly circular or irregular (few edge-on disks). For a sanity check, we also run GALFIT (Peng et al. 2002) for galaxies in the sample. We fit our galaxies with the Sersic profiles and a free Sersic index. Giving the Source Extractor parameters as an initial guess to GALFIT, we get half-light radii that are consistent with those from Source Extractor. We compare the half-light radii derived from Source Extractor and GALFIT for a subsample of our galaxies in Figure 1. The errors derived from the GALFIT algorithm are very small, with a median value of 2%.

Zoom In Zoom Out Reset image size

To be consistent with El-Badry et al. (2016), we only keep galaxies with axis ratios > 0.5 and stellar masses between 10⁷ and 10^9.6 M_⊙, where burstiness is predicted to be more prominent compared to other mass ranges (El-Badry et al. 2016). This leaves us with a final sample of 86 galaxies.

We present our size measurements as a function of stellar mass in Figure 2. The logarithm of half-light radius (R_e ) in our LVL sample linearly increases with the logarithm of mass with a slope of 0.282 ± 0.025. Our R_e values are in agreement with the r-band effective radii of SDSS z ∼ 0 galaxies measured in the NASA-Sloan Atlas (Blanton et al. 2011).

Zoom In Zoom Out Reset image size

At fixed mass, there is ∼0.8 dex scatter inferred by the standard deviation. This scatter is the main driver of this work and is predicted to be partially driven by the feedback-driven outflows generated by episodic bursts of SF (El-Badry et al. 2016). Here we aim to examine this prediction and see whether there is any relationship between the size fluctuation and SF change in the observed low-mass galaxies.

3.1. Observed Size–Luminosity Relation

In order to test the relation between the size and the SFR (or, equivalently, the luminosity) of the sample, we need to first subtract off the mass dependency from the R_e distribution. In particular, at any given mass, we are interested in each galaxy's deviation in R_e from the average trend shown in Figure 2. For this, we introduce a new definition of size, which is the ratio of R_e to the average R_e (${R}_{e}/\overline{{R}_{e}}$). Technically, at each point we divide the observed R_e by the average $\overline{{R}_{e}}$ derived from the line fit through the log(R_e )–log(M_*) relation and refer to it as normalized-to-average size.

Figure 3 shows postage stamps of R-band images of a sample of LVL galaxies that are scaled by their distances as if they are all at the same distance and then are sorted by their ${R}_{e}/\overline{{R}_{e}}$. The top left panel shows the most compact galaxy in the sample with the smallest ${R}_{e}/\overline{{R}_{e}}$, while the bottom right panel indicates the galaxy with the largest ${R}_{e}/\overline{{R}_{e}}$.

Zoom In Zoom Out Reset image size

Next, we calculate the dust-corrected Hα and far-UV (FUV) luminosities of the sample (Weisz et al. 2012). Hα and FUV luminosities are the two common probes of the SFRs in galaxies with timescales of 5 and 10–100 Myr respectively, depending on the SF systems (Kennicutt 1998; Emami et al. 2019; Flores Velázquez et al. 2021).

In order to take out the mass dependency from the Hα (and FUV) luminosity distributions, we divide the Hα (FUV) luminosity at a given mass by the average luminosity derived from the line fit through the log L_Hα (log L_UV) and log M_* relation at that mass, the same procedure we used to derive the normalized-to-average size values (${R}_{e}/\overline{{R}_{e}}$). We then split the sample into four mass bins with equal number of galaxies in each bin and investigate the size–luminosity relation in each mass bin assuming that galaxies of similar stellar masses likely share common SF properties (Weisz et al. 2012; Emami et al. 2019). In particular, lower-mass galaxies should have a lower potential well and thus may be more likely to show signs of episodic SF in their sizes.

Figure 4 shows the logarithm of normalized-to-average Hα and FUV luminosities, namely, L_Hα /$\overline{{L}_{{\rm{H}}\alpha }}$ and L_FUV/$\overline{{L}_{\mathrm{FUV}}}$ as a function of log(${R}_{e}/\overline{{R}_{e}}$) for four mass bins. We find a positive correlation between log(L_FUV/$\overline{{L}_{\mathrm{FUV}}}$) and log (${R}_{e}/\overline{{R}_{e}}$) at masses < 10^7.85 M_⊙. The correlation becomes weaker and more scattered toward larger mass bins. For the log(L_Hα /$\overline{{L}_{{\rm{H}}\alpha }}$)–log(${R}_{e}/\overline{{R}_{e}}$) relation, we also see a positive slope but less steep and more scattered than the log(L_FUV/$\overline{{L}_{\mathrm{FUV}}}$)–log (${R}_{e}/\overline{{R}_{e}}$) relation. We report the slopes of the best-fit relationships for each mass bin and the Spearman rank correlation coefficients in Table 1. A larger Spearman rank correlation coefficient indicates a stronger correlation between two quantities.

Zoom In Zoom Out Reset image size

Table 1. The Slopes and Spearman Rank Correlations (ρ) of the Logarithm of FUV and Hα Normalized-to-average Luminosities versus Logarithm of Normalized-to-average Half-light Radius (${R}_{e}/\overline{{R}_{e}}$) for Four Stellar Mass Bins in the LVL Sample

$\mathrm{log}({M}_{* }/{M}_{\odot }$)	log(LFUV/$\overline{{L}_{\mathrm{FUV}}}$)–log(${R}_{e}/\overline{{R}_{e}}$)	log(LHα /$\overline{{L}_{{\rm{H}}\alpha }}$)–log(${R}_{e}/\overline{{R}_{e}}$)
	slope (ρ)	slope (ρ)
7–7.5	1.00+/−0.28 (0.66)	0.70+/−0.40 (0.37)
7.5–7.85	0.91+/−0.31 (0.65)	0.50+/−0.50 (0.37)
7.85–8.5	0.26+/−0.32 (0.03)	0.03+/−0.33 (−0.06)
8.5–9.6	0.33+/−0.34 (0.24)	0.09+/−0.32 (0.05)

Download table as:  ASCIITypeset image

In order to better understand these observed trends, we compare to dwarf galaxies in hydrodynamical zoom-in simulations.

In Section 3, we found that for masses below 10^7.85 M_⊙ the normalized-to-average size is correlated with the FUV normalized-to-average luminosity at fixed stellar mass. This supports predictions from hydrodynamical simulations (El-Badry et al. 2016) suggesting that the size fluctuation and stellar migration in low-mass systems are consequences of feedback-driven outflows, which in turn are driven by the cumulative effects of many starburst episodes over time.

Based on this argument, SF is triggered when cold gas accretes into the center of the galaxy and forms stars. Stellar feedback then creates a galactic wind and drives gas to larger radii, resulting in a reduction in SF. Since the velocity of the expelled gas does not exceed the escape velocity (Muratov et al. 2015), it soon cools, re-accretes back to the center, and resumes SF. This periodic gas displacement across the galaxy creates fluctuations in the central gravitational potential, which in turn leads to variation in the stellar distribution and the measured size of the galaxy. This theoretical picture can explain the relationship between the FUV luminosity and size in our observed low-mass sample.

However, it is not yet understood why this dependency is less effective between the size fluctuation and short-timescale SFR indicators, i.e., Hα luminosity. Since Hα and FUV luminosities trace the SFR on different timescales (Flores Velázquez et al. 2021), a detailed investigation on the timescale over which the galaxy size fluctuation and SF variation occur is critical. In this section we aim to compare the SFHs with the half-light radius time sequences for a set of simulated galaxies in order to find quantitative explanations for our observed findings.

Feedback in Realistic Environments (FIRE-1, FIRE-2; Hopkins et al. 2014, 2018) is a set of simulations with models for SF and stellar feedback implemented within cosmological baryonic simulations. In the FIRE model, energy and momentum inputs from stars are taken from stellar evolution calculations and models of the unresolved Sedov–Taylor phase, without tunable parameters. The simulations have been successful in reproducing several key observables, including realistic galactic outflows (Muratov et al. 2015, 2017), the dense H i content of galaxy halos (Faucher-Giguère et al. 2015), the mass–metallicity relation (Ma et al. 2016), the mass–size relation and age/metallicity gradients (El-Badry et al. 2016), cored dark matter profiles (Chan et al. 2015; Oñorbe et al. 2015), stellar kinematics (Wheeler et al. 2017), the Kennicutt–Schmidt relation (Orr et al. 2018), observed abundance distributions (Escala et al. 2018), and a realistic population of satellites around Milky Way–mass hosts (Wetzel et al. 2016).

Here we use FIRE-2 simulations of eight isolated low-mass galaxies spanning over a mass range of 10⁷–10^9.6 M_⊙, analogous to our local observed sample, and explore their SF and size evolution over the past 2 Gyr (z ∼ 0.15–0).

4.1. Simulated Size–SFR Relation

To calculate the R-band half-light radii of the simulated galaxies, we need to first generate mock R-band images of each snapshot in the simulation. The time spacing between consecutive snapshots is ∼20 Myr. Following procedures in Chan et al. (2018), we generate a table of R-band luminosity-to-mass ratios of single stellar populations for a range of metallicity and age. For that we run the Flexible Stellar Population Synthesis (Conroy et al. 2009, 2010) model and assume the latest Padova stellar evolutionary tracks (Marigo & Girardi 2007; Marigo et al. 2008), Kroupa initial mass function (Kroupa 2001), and MILES stellar library. We determine the R-band luminosity-to-mass ratio of each simulated star particle by interpolating the generated table at the particle's age and metallicity and then multiply the resulting luminosity-to-mass ratio by the mass of the star particle to get the R-band luminosity. We then project the R-band luminosities of each simulation output over a 20 kpc × 20 kpc region onto 1400² mesh grids and generate mock R-band images. We have examined our analysis on a larger number of mesh grids and found our size measurements robust to spatial resolutions. Since most of these galaxies are spherical, the size measurements are independent from the orientation we choose to view the galaxy. We show snapshots of R-band mock images of galaxy m11c in Figure 5, which indicate that the galaxy undergoes a significant size fluctuation across ∼150 Myr.

Zoom In Zoom Out Reset image size

We calculate the half-light radius for each snapshot in all simulated galaxies by passing the R-band mock images to Source Extractor and use the same input parameters as those used for our observed galaxies (See Section 2). We also use a point-spread function (PSF) equivalent to the median seeing of our observed sample (14) assuming that the simulated galaxies are at a 5 Mpc distance (which is the median distance of the observed sample) and find that the chosen PSF is much smaller than the galaxy sizes. In order to account for the potential impact of sky noise on our size measurements, we add stochastic sky backgrounds with an average surface brightness ∼26 mag arcsec⁻² to our mock images. We find that the differences in the estimated half-light radii with and without sky backgrounds are small and do not affect our analysis. Figure 6 shows the R-band half-light radius as a function of stellar mass for FIRE-2 galaxies. The distribution of the simulated galaxies is consistent with that of the observed LVL galaxies in the R_e –mass parameter space. Also, similar to our observed sample, we find our size measurements to be in agreement with those of z ∼ 0 SDSS galaxies (Blanton et al. 2011). We also see that R_e fluctuates by more than a factor of 2 over ∼200 Myr in most of the galaxies.

Zoom In Zoom Out Reset image size

4.2. Size Fluctuation as a Consequence of Bursty Star Formation

Next, we compare the SFHs with the R_e time sequences of the simulated galaxies. Figure 7 shows the SFH and R_e time sequences for galaxy m11q (10^8.6 M_⊙). The local minima and maxima in the R_e time sequences coincide with certain phases in the SFH; R_e is at its lowest when SF has just begun rising (nearly after a few Myr; green arrows), and it peaks once SF shuts down (purple arrows). The physical interpretation is that when SF has just begun rising, gas in the center is highly concentrated and forms stars at small radii, and that the majority of the stellar mass is in the center. After supernovae begin to explode, the energetic outflowing winds drive away the gas, resulting in a shallow central gravitational potential that drives stars to outer radii. The subsequent supernovae cumulatively drive gas and stars farther away, and the galaxy expands and reaches its maximum size. Once SF has completely shut off, the central potential returns to a relaxation phase, the cooled gas and stars reaccrete, and the galaxy starts to contract until the next burst cycle begins. We see this phenomenon in all of our simulated galaxies at a mass range of 10⁷–10^9.6 M_⊙. Based on this argument and Figure 7, it can be inferred that the size of the galaxy traces the SF changes with a time delay.

Zoom In Zoom Out Reset image size

4.3. What Is the Time Delay from SF Change to Galaxy Size Change?

In order to determine how long it takes the size of a galaxy to respond to the SF change, i.e., the time delay from SF change to the galaxy size change, we calculate the cross-correlation between the R_e time sequence and SFH for each galaxy in the FIRE-2 sample. Since our R_e data points are spaced by 20 Myr while our SFRs are sampled at every single Myr, it is very important to first smooth the time sequence of the latter to a time resolution close to that of R_e and then sample the time sequences for 20 Myr intervals; otherwise, it is likely that the long-term SF fluctuations will be lost behind the instantaneous SF fluctuations.

In the cross-correlation algorithm, the correlation between R_e and SFR time sequences with lag k is defined as

$$\begin{eqnarray}&&\sum _{n}{R}_{e}[n+k].{\mathrm{SFR}}^{* }[n],\end{eqnarray} \tag{ 1 }$$

where SFR^*[n] is the complex conjugate of SFR at the nth snapshot and R_e [n + k] is the R_e value at the (n+k)th snapshot. Figure 8 shows the cross-correlation between R_e and SFR as a function of time lag for galaxy m11q. The peak in the cross-correlation plot indicates the time lag between the SFH and R_e time sequences, by which the R_e time sequences correlate best with the SFH if the SFH is shifted by that amount. We show the same plots for the rest of the sample in Figure 11 in the Appendix. The time lag at the maximum cross-correlation is reported in Table 2. In six out of eight galaxies, the cross-correlation peaks at ∼50 Myr time lag. This implies that it takes the galaxy's size and the stellar migrations roughly 50 Myr to react to any changes in the SFR. For the other two galaxies (m11b, m11h), the cross-correlation peaks at zero Myr time lag between R_e and SFR. This is because in these two galaxies the starbursts occur immediately one after another such that there is no quenching time between consecutive bursts, making it difficult for the cross-correlation algorithm to identify individual extrema in the SFH and R_e time sequences.

Zoom In Zoom Out Reset image size

Table 2. Time Lag at the Maximum Cross-correlation between SFR and Half-light Radius (R_e ) in FIRE-2 Galaxies

Name	log(M*/M⊙)	Re –SFR Time Lag (Myr)
	z ∼0	(Cross-correlation)
m10y	7.03	48
m11a	8.0	52
m11b	8.0	0
m11q	8.6	50
m11c	8.9	63
m11i	8.95	43
m11h	9.5	0
m11d	9.56	56

Download table as:  ASCIITypeset image

We also perform the cross-correlation algorithm on three SFR observables with different timescale sensitivities, including the Hα luminosity and FUV (1500 Å) and near-UV (NUV; 3500 Å) luminosity densities. Our goal is to see which one of these SFR observables has the shortest time lag with the R_e and thus is better (positively) correlated with the galaxy size. The medians of the time lags between R_e and Hα, FUV, and NUV for all the galaxies are 42, 28, and 8 Myr, respectively. This means that R_e and NUV luminosity are more correlated with each other (compared to the R_e with Hα or FUV) and NUV luminosity traces the SFR on a comparable timescale to the R_e changes. This is expected given that the FUV and Hα light comes from short-lived, massive, and hot stars that trace the SF change on a short timescale, while NUV light comes from longer-lived, less massive stars that take a longer time to trace the SF variations.

In Figure 9 we show the normalized R_e time sequences along with the normalized SFH that lags by 50 Myr (SFR₅₀), as well as the normalized NUV, FUV, and Hα time sequences for galaxy m11q. It is evident that an SFH with a 50 Myr time lag is well in phase with the R_e time sequence. The correlation between R_e and the three SFR observables decreases from NUV, to FUV, to Hα luminosity, respectively, until the R_e and Hα time sequences become fully uncorrelated. This can also be seen in almost all of the FIRE-2 galaxies (see online journal for the rest of the FIRE-2 sample).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (7 panels)

Tarball of figure set images

In order to directly compare these results from simulated galaxies with our observed results in Section 3, we plot all the SFR indicators as a function of R_e in Figure 10 and fit a line to each distribution. We present the slopes and the Spearman rank correlation coefficients (ρ) in Table 3. It is evident that there is a positive slope in the R_e –SFR₅₀ relation and the slope decreases from left to right, ultimately reaching an uncorrelated R_e –L_Hα . We also see a shallow FUV–R_e correlation in most of the galaxies, which is in agreement with our observed findings in Figure 4. Unlike in the observations, however, we do see a correlation at all masses. This is in tension with what we found in our observed sample, in which the correlation was stronger at masses below 10^7.85 M_⊙ and became less correlated toward higher masses. This is due to the bursty nature of SF in FIRE simulations, in which the timescale of SF bursts is much shorter than the timescale implied by observations of real galaxies at high mass (Muratov et al. 2015; Sparre et al. 2017; Emami et al. 2019; Stern et al. 2021). Furthermore, we see that in Figure 10, for all four SFR indicators (SFR₅₀, L_NUV, L_FUV, and L_Hα ), the slopes are biased by a few points that show very low SFR. These points are associated with the post-burst phases that are followed by long periods of no SF. In these cases, SFR decreases rapidly, while R_e takes a longer time to decrease and reach an equilibrium. One example can be seen in Figure 9, for which during the quenching time at 13.55 Gyr and the subsequent burst at 13.75 Gyr SFR decreases rapidly while R_e decreases at a much lower rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Slopes and the Spearman Rank Correlation (ρ) of the log(R_e )–log(SFR) Relation for Different SFR Indicators in Each FIRE-2 Galaxy

Name	Re –SFR50	Re –LNUV	Re –LFUV	Re –LHα
	Slope (ρ)	Slope (ρ)	Slope (ρ)	Slope (ρ)
m10y	1.28 ± 1 (0.03)	0.62 ± 1.38 (−0.08)	0.12 ± 1.56 (−0.11)	−4.22 ± 2.25 (−0.30)
m11a	2.61 ± 1.05 (0.32)	0.59 ± 0.47 (0.12)	−0.07 ± 0.86 (−0.04)	−3.47 ± 1.73 (−0.28)
m11b	0.70 ± 0.24 (0.33)	0.59 ± 0.19 (0.32)	0.41 ± 0.27 (0.17)	0.37 ± 0.46 (0.12)
m11q	3.94 ± 1.16 (0.67)	0.88 ± 0.44 (0.38)	0.78 ± 0.80 (0.24)	−1.80 ± 1.73 (−0.08)
m11c	3.20 ± 1.22 (0.39)	2.53 ± 0.55 (0.48)	1.80 ± 0.76 (0.17)	−1.79 ± 1.51 (−0.16)
m11i	1.93 ± 0.44 (0.44)	0.77 ± 0.23 (0.33)	0.72 ± 0.38 (0.16)	0.11 ± 0.61 (−0.04)
m11h	0.93 ± 0.29 (0.31)	1.12 ± 0.25 (0.45)	0.84 ± 0.27 (0.28)	0.81 ± 0.38 (0.21)
m11d	2.54 ± 0.66 (0.44)	1.06 ± 0.25 (0.38)	1.59 ± 0.47 (0.29)	0.30 ± 1.00 (0.03)

Download table as:  ASCIITypeset image

5.1. Previous studies

5.2. The Role of Galaxy Formation Histories and Mergers in the ${R}_{e}/\overline{{R}_{e}}$–SFR Relation

5.3. Implications for High-redshift Tests with the James Webb Space Telescope

We explored the effect of bursty SF on galaxy sizes, which results from the implementation of stellar feedback onto the dark-matter-only simulations. Winds driven by outflows induce fluctuations on the gravitational potential and cause a radial stellar migration in low-mass galaxies. We tested this prediction on samples of both observed local dwarf galaxies and simulated late-time dwarfs across a mass range of 10⁷–10^9.6 M_⊙, where theoretical simulations suggest feedback to be most efficiently altering galaxies' kinematics (El-Badry et al. 2016). For our observed sample, we used Hα and FUV luminosities (L_FUV) to probe different timescales of the SFR and the R-band images to measure the half-light radii (R_e ). We found that the half-light radius correlates well with L_FUV while the correlation becomes weaker between R_e and L_Hα . Furthermore, the R_e –L_FUV relation is much more significant at masses below ∼ 10⁸ M_⊙ and becomes shallower at higher masses. These findings suggest that the half-light radius changes owing to the SF variations and reacts to this variation within a timescale similar to that of FUV luminosity (i.e., 10 Myr or longer), but not close to the Hα timescale at all, which is a few Myr. To understand the underlying physics of this size–SFR correlation and determine the approximate timescale over which the size traces the SFR variations, we explored eight low-mass galaxies from the FIRE-2 cosmological simulations at late cosmic times (over the past 2 Gyr) that mimic our observed local sample. We summarize our findings as follows:

• 1.  
  By visual inspections of the half-light radius time sequence and the SFH of galaxies, it is evident that there is a delay from the onset of a burst of SF to when half-light radius increases. The physical interpretation is that when SF begins rising, gas and stars are highly concentrated and the galaxy is compact. After subsequent supernovae, the resulting outflows drive gas and stars to larger radii and the galaxy expands. This expansion continues until the SF completely shuts down and then cooled gas and stars migrate back to the center.
• 2.  
  In order to assess the timescale by which R_e traces the SF change, we performed the cross-correlation technique. We found that there is a 50 Myr time lag between the SFH and half-light radius in almost all galaxies, suggesting that the size of the galaxy traces the SF changes with a 50 Myr delay. This is also evident as a positive slope in the half-light radius–SFR₅₀ relation for almost all the FIRE galaxies.
• 3.  
  Using multiple common SFR observables, including the Hα, FUV, and NUV luminosities, and cross-correlating these with the half-light radius, we found that the time sequence of NUV luminosity is more correlated with the half-light radius time sequence, and the correlation decreases for the FUV and then Hα time sequences, respectively. This is consistent with our observed results, for which half-light radius is better correlated with the FUV luminosity than the Hα luminosity.
• 4.  
  The half-light radius–SFR correlation exists in all the FIRE galaxies spanning a mass range of 10⁷–10^9.6 M_⊙. This is in tension with what we found in our observed sample, in which the correlation can only be seen in masses below 10^7.85 M_⊙. This is likely related to the bursty nature of most massive galaxies in the FIRE simulations at low redshifts (Sparre et al. 2017; Emami et al. 2019).
• 5.  
  With the high surface brightness sensitivity of the near-IR imaging cameras on the James Webb Space Telescope, we will be able to test the effect of bursty SF on the size of the high-redshift galaxies, while the biases against low surface brightness galaxies are significantly minimized.

We thank the anonymous referee for providing constructive comments that helped improve the quality of this paper. We also thank T. K. Chan and Amanda Pagul for their helpful discussions. This material is based on work supported by UC Riverside Overhead Return Grant under grant No. A01868-19900-40-FPBS.

Figure 11 shows the cross-correlation between Re and SFR as a function of time lag for FIRE-2 galaxies same as the plot shown in Figure 8 for galaxy m11q. The peak in the cross-correlation plots indicate the time lag between the SFH and Re time sequences, by which the Re time sequences correlate best with the SFH if the SFH is shifted by that amount. The time lag at the maximum cross-correlation is reported in Table 2. In six out of eight galaxies, the crosscorrelation peaks at ∼ 50 Myr time lag. This implies that it takes the galaxy's size and the stellar migrations roughly 50 Myr to react to any changes in the SFR.