For the analysis of the ionized gas morphology and kinematics, we use NIRCam grism data in the F444W filter (m) obtained with FRESCO (P. A. Oesch et al. 2023) in both GOODS fields and in the F356W filter (m) obtained with CONGRESS (PIs: Egami, Sun, PID: 3577) in GOODS-N. These filters are chosen as they probe the emission () at . The observations for both surveys are taken using the row-direction grisms (GRISMR mode) with modules A and B, with a spectral resolution of across the filters used. The reduction technique used for the grism data is detailed in F. Sun et al. (2023) and J. M. Helton et al. (2024). The continuum-subtracted emission-line maps used for the morpho-kinematic modelling (Section 2.6) are obtained by applying row-by-row median filtering using a 2-step iterative technique as described in D. Kashino et al. (2023).

The FRESCO and CONGRESS grism programmes also obtained imaging data in the F182M and F210M, and F090W and F115W filters, respectively. In addition, there is a wealth of imaging data in the GOODS fields from the JADES survey (D. J. Eisenstein et al. 2026; M. J. Rieke et al. 2023b), a Guaranteed Time Observations (GTO) programme of the NIRCam and NIRSpec instrument teams that obtained imaging of an area of arcmin with 8–10 filters. This overlap provides additional imaging in the wide bands F090W, F115W, F150W, F200W, F277W, F356W, and deeper F444W observations, and medium bands F335M and F410M observations. For certain regions in our sample, we also have additional medium bands F182M, F210M, F430M, F460M, and F480M from the JWST Extragalactic Medium Band Survey (JEMS; C. C. Williams et al. 2023) in GOODS-S. This wide range of filters is not only optimal to probe sizes across different wavelengths, but also to provide tighter constraints on the SED modelling (Section 2.4). We use high-resolution data and photometric catalogues obtained from drizzled mosaics, and the full details on the data reduction, catalogue generation, and photometry can be found in M. J. Rieke et al. (2023b); B. Robertson et al. (2024).

Our sample is taken from A. L. Danhaive et al. (2025), where a full description of the sample selection can be found, but we offer a summary here. The sample of emitters from the FRESCO and CONGRESS surveys is built using photometric redshift estimates from eazy (G. B. Brammer, P. G. van Dokkum & P. Coppi 2008; K. N. Hainline et al. 2024), which are then verified through Gaussian modelling of the 1D integrated spectrum containing the line, amongst others. The derived spectroscopic redshifts are then confirmed through visual inspection. This process resulted in a catalogue of   emitters with S/N > 5 detections, which is in good agreement with previous work using the same grism data (A. Covelo-Paz et al. 2025). In this work, we make a more conservative cut of , and obtain 393 galaxies in CONGRESS and 189 in FRESCO (J. M. Helton et al. 2024; X. Lin et al. 2026).

However, in order to recover reliable constraints on the kinematics from grism data, we need not only high enough S/N, but also a favourable galaxy orientation on the sky. In fact, if the major axis is aligned to the grism dispersion direction, the effects of the velocity and velocity dispersion are completely degenerate with each other, and with the size and morphological profile of the galaxy. We therefore require that the position angle (PA) of the galaxy, defined from the positive y-axis, be distinct from the dispersion direction (defined by ). To achieve this, we only select galaxies with . We also remove objects that have multiple components and therefore are not well modelled by a single Sérsic profile, such as merging systems. The PA and S/N cuts can be more or less conservative, resulting in the three subsamples presented in A. L. Danhaive et al. (2025):

In this paper, we use the full sample of 213 galaxies in our analyses and make no further distinction between the three subsamples. This distinction was important in A. L. Danhaive et al. (2025) where the focus is on kinematics, but in this morphology-centred paper, we prioritize sample size over tight kinematic constraints.

As detailed in A. L. Danhaive et al. (2025), we infer stellar masses, SFRs averaged over different time-scales, and constraints on the dust attenuation from SED modelling with Prospector (B. D. Johnson et al. 2021). We fit all of the available photometry from JADES (Section 2.2), including MIRI imaging in F770W and F1280W, in conjunction with additional bands from JEMS and FRESCO (F182M, F210M, F430M, F460M, and F480M), and available photometry from HST’s Advanced Camera for Surveys (ACS) bands F435W, F606W, F775W, F814W, and F850LP (the GOODS survey M. Dickinson et al. 2004; M. Giavalisco et al. 2004) released through the Hubble Legacy Field programme (K. E. Whitaker et al. 2019). Following S. Tacchella et al. (2023), we model the photometry self-consistently with the available line fluxes from the FRESCO and CONGRESS surveys, and, for FRESCO galaxies in the GOODS-N field, we also include, when detected, Å, Å, H, and [S iii]9533 Å fluxes from the CONGRESS survey. The redshift in the SED fitting is fixed to the spectroscopic redshifts derived from the grism spectra, which allows us to break the degeneracy between redshift and physical properties such as stellar mass.

We present our sample in Fig. 1, where we show the stellar mass as a function of redshift. We colour-code our sample by sSFR, which is defined as  and is a good probe of where galaxies lie with respect to the star-forming main sequence (SFMS). It is apparent that, at low masses , we are biased to highly star-forming systems with high sSFRs. This is a direct result of our integrated flux S/N cut. The sSFR measures the SFR per unit stellar mass of a galaxy, and it roughly equals the inverse of the mass doubling time-scale . This implies that the galaxies in our sample, if they keep a constant sSFR, would double their stellar masses in only . Although galaxies oscillate around the SFMS and do not have constant sSFRs, this still highlights the short time-scales on which galaxies grow at .

In order to constrain the sizes of the galaxies in our sample at different rest-frame wavelengths, we use the array of available wide-band imaging from the JADES survey. For galaxies at (), the rest-frame FUV ( Å), NUV ( Å), and optical ( Å) continuum are well-probed by the F090W, F150W, and F277W (F090W, F200W, and F356W) filters (Table 1). The imaging in all of these filters has a comparable exposure time ( 5h) for most of the galaxies in our sample. The sizes are constrained using morpho-kinematic modelling of the F356W and F444W grism data, respectively (see Section 2.6).

We fit the relevant wide band images with the fully Bayesian code pysersic (I. Pasha & T. B. Miller 2023). pysersic models imaging data from any filter to infer the best-fitting Sérsic profile(s) parameters. pysersic forward-models the image and accounts for the PSF convolution, so that the resulting morphological parameters are intrinsic, meaning PSF-deconvolved. We approximate the PSF in each filter using the model PSFs (mPSFs) constructed by mosaicking WebbPSF models repeatedly over the field identically to our exposure mosaics and then measuring the average PSF. This method is presented in Z. Ji et al. (2026). The model PSFs are then resampled down to the specific filter pixel resolution and used for modelling. In terms of sampling, we use the stochastic variational inference (SVI)-flow mode of pysersic, where the posterior distributions of the model parameters are inferred using normalizing flows. We note that this method is consistent with results for the same objects obtained using the pysersic Markov Chain Monte Carlo (MCMC) sampling mode (C. Carreira et al. 2026).

We model our galaxies with a one-component Sérsic profile (J. L. Sersic 1968):

and derive measurements of position angle , ellipticity e, Sérsic index n, light centroid , and, importantly, half-light radius . Each band is independently fit. Additionally, we fit for a flat sky background and we mask any contaminant source from the cut-out. While masking neighbouring sources is common in morphological modelling (e.g. T. B. Miller et al. 2025; N. Allen et al. 2025), another established method replaces masking with explicit modelling of any additional sources within the cut-out (e.g. K. Ormerod et al. 2024; H. M. O. Stephenson et al. 2025). The derived best-fitting parameters can change depending on which method is used, but they remain consistent within .

All of the free model parameters have flat priors, whose bounds are purely related to physical constraints. For instance, the total amplitude and the size need to be positive, the ellipticity has to be between 0 and 1, and the Sersic indices can span the full range between 0.36 and 8. Fig. 2 shows the fits for a galaxy in our sample, as an example for our morphological modelling procedure. Despite some underlying structure present in the residuals, we obtain a satisfactory fit given our simple model. Specifically, this structure is not expected to strongly affect the measurement of sizes, which is dominated by the main bright component at high-redshift (e.g. M. Mosleh, R. J. Williams & M. Franx 2013).

Using the geko code, we forward-model a specific emission line, in this case , morphology and kinematics and fit the NIRCam grism data using a Bayesian inference framework. The full methodology is detailed in A. L. Danhaive et al. (2025) and the full code is presented in A. L. Danhaive & S. Tacchella (2025),¹ but we summarize the key elements here.

The emission line morphology is modelled with a single Sérsic profile with the six free parameters outlined in Section 2.5. To generate a prior for these parameters, needed to help alleviate the morphology-kinematics degeneracy of the 2D grism data, we use results from our pysersic modelling of the F150W band, which probes the UV continuum in both redshift ranges. For galaxies not covered by the F150W filter, we use the nearby F182M filter. The width of the priors is estimated by doubling the uncertainties derived with pysersic to account for the fact that we are not measuring the true emission-line morphology, but stellar continuum emission which probes similar star formation time-scales as . Importantly, the tests presented in A. L. Danhaive & S. Tacchella (2025) show that the recovery of sizes is robust even when the morphological priors are offset from their true values, as well as for a wide range of different kinematic configurations. Conversely, the recovery of the true kinematic parameters is successful for most morphological configurations, except for the most compact cases.

The morphological model is then convolved with the kinematics before being convolved with the instrument’s PSF and line-spread function (LSF) and projected in the grism 2D space. The kinematics are modelled with an arctangent velocity curve (S. Courteau 1997; S. H. Miller et al. 2011)

and a constant velocity dispersion. The inclination of the velocity curve is set to the value derived from the morphological fit. Specifically, because galaxies are thought to be thicker at high-redshift (A. der Wel et al. 2014b), we compute the inclination assuming an intrinsic disc axis ratio, which is the ratio of scale height to scale length, (S. Wuyts et al. 2016; R. Genzel et al. 2017; S. H. Price et al. 2020; H. Übler et al. 2024), where corresponds to the thin disc approximation. The free kinematic parameters are hence: , the asymptotic velocity, , the turn-around radius, , the position angle which is fit independently from the morphological one, and intrinsic velocity dispersion . The kinematic parameters have uniform priors, with forced to be below , motivated by lower redshift studies finding (S. H. Miller et al. 2011). The prior for the kinematic position angle is the same as for , despite being fit independently.

The SMR is typically fit by a single power law in the literature (e.g. A. der Wel et al. 2014a; L. A. Mowla et al. 2019). Following these works, we adopt a linear relation in logarithmic space, which accounts for both the redshift and the mass dependence of galaxy sizes:

where is the half-light radius, is the stellar mass, and is the slope of the mass dependence, is the normalization (the size for a galaxy at ), and parametrizes the redshift dependence. We favour this method over fitting a redshift-independent equation in different redshift bins (e.g. N. Allen et al. 2025; L. Yang et al. 2025) since our sample spans a relatively narrow redshift range, corresponding to a cosmic time interval of Gyr, and suffers from a decrease in galaxies at higher redshifts (see Fig. 1). Furthermore, fitting for both the mass and the redshift dependence simultaneously allows us to study the effect of each factor independently. This is difficult to do otherwise, since the average mass of galaxies naturally decreases with redshift. We also fit for the intrinsic scatter around the SMR, assuming that the sizes follow a lognormal distribution at fixed mass and redshift, .

In order to infer the best-fitting parameters of the SMR, and their uncertainties, we use the Bayesian inference package emcee (D. Foreman-Mackey et al. 2013). emcee samples the parameter space using an implementation of the affine invariant ensemble sampler for MCMC (J. Goodman & J. Weare 2010), resulting in posterior distributions for the model parameters. We adopt uniform priors for all parameters. The uncertainties in this work are defined with the and quantiles of the posterior distributions.

We present our results for the , FUV, NUV, and optical SMR (equation 4) at in Fig. 3, using our sample of 213 emitters from the FRESCO and CONGRESS surveys. The sizes are derived from the morpho-kinematic fitting of the grism data (Section 2.6), and the stellar continuum sizes are obtained from modelling with pysersic (Section 2.5). We first fit the SMR only for galaxies with , where our sample is most representative in terms of mass completeness (). In fact, our selection bias comes from the integrated S/N limit, but since luminosity scales with stellar mass at fixed sSFR, this is approximately equivalent to a stellar mass cut. Above , we find a flat trend between sSFR and stellar mass, suggesting completeness, whereas below we are biased to high sSFR systems. About half of our sample (58 per cent) is represented in the mass cut.

In order to quantify the effects of this choice, we also include a fit done using the full sample. The best-fitting parameters for the ‘fiducial’ SMR and the one derived from the full sample are summarized in Table 2. We show the posterior distributions for the fiducial , FUV, NUV, and optical fits in Fig. A1. We plot both best-fitting SMRs at on Fig. 3, where we also compare to other SMRs at similar redshifts from the literature (T. B. Miller et al. 2025; N. Allen et al. 2025; L. Yang et al. 2025). We also compare our observed SMR with predictions from simulations, namely thesan (E. Garaldi et al. 2022; R. Kannan et al. 2022; A. Smith et al. 2022), presented in X. Shen et al. (2024), thesan-zoom (R. Kannan et al. 2025), presented in W. McClymont et al. (2025b), tng50 (D. Nelson et al. 2019; A. Pillepich et al. 2019), presented in L. Costantin et al. (2023), and flares (C. C. Lovell et al. 2021; A. P. Vijayan et al. 2021), presented in W. J. Roper et al. (2022).

We find that the best-fitting SMRs with and without the mass cut are consistent within the uncertainties for the FUV, NUV, and optical sizes, although the full sample slopes are all shallower by dex. This effect is stronger for the SMR which is almost flat when the full sample is included (). The posterior distributions for the full sample fit are shown in Fig. A2. This flattening could be partially caused by the large scatter of Sérsic indices found at low masses, where we find preferentially higher values of n (Fig. A4). Because of the degeneracy between n and (equation 2), higher Sérsic indices naturally result in larger sizes. Also, as highlighted in Fig. 1, most of the low-mass galaxies in our sample lie above the main sequence, , where

and is the SFR, averaged over the last 10 Myr, for a galaxy of mass on the main sequence. We adopt the SFMS from C. Simmonds et al. (2025). As we will discuss in Section 5, this bias could cause the observed flattening of sizes if galaxies above the main sequence have boosted sizes. In the following sections, we focus on our fiducial fit, unless otherwise specified, but note that the main conclusions remain consistent with the full sample fit.

We begin by focusing on the normalization of the SMR, parametrized by , which we measure at and (equation 4). The sizes, which probe the most recent star formation (), are the largest, with kpc at . At this stellar mass, the stellar continuum sizes are all below 1 kpc, with kpc and kpc. At all wavelengths, we see a large scatter of sizes at fixed stellar mass, spanning dex. The larger sizes are typically associated with the inside–out growth of galaxies, which we will discuss further in Section 4, as they suggest that the recent star formation extends further out than the older stellar populations probed by the rest-frame UV and optical. Interestingly, we do not find a significant difference between UV and optical sizes, consistent with other studies from JWST (T. Treu et al. 2023; T. B. Miller et al. 2025; N. Allen et al. 2025). This could be due to a number of reasons. In terms of observational effects, outshining by young stellar populations could cause all of the continuum sizes to be probing stars of similar (young) ages. Also, dust attenuation can play a role, as we discuss in Section 4.3. Physically, the lack of difference in continuum sizes could imply that young stars are not preferentially forming in the outskirts of galaxies.

In Fig. 4, we put our measured galaxy sizes in the context of other works from the literature. Compared to J. Matharu et al. (2024), who compute sizes at similar redshifts based on stacks of emitters in FRESCO, we find larger sizes by dex. Stacking suppresses features that are not perfectly aligned, so both small misalignments and structural variations could contribute to this discrepancy, particularly given that the maps are not corrected for kinematic distortions. Finally, although J. Matharu et al. (2024) account for PSF effects in their morphological modelling, they do not take into account the grism LSF for the stacks, as we do in this work, which could also affect their size measurements. In H. M. O. Stephenson et al. (2025), the emission is traced by narrow-band imaging for 23 star-forming galaxies at as part of the JWST Emission Line Survey (JELS; K. J. Duncan et al. 2025; C. A. Pirie et al. 2025). As shown on Fig. 3, the median sizes for their sample are dex smaller than our fit at . This discrepancy could be driven by biases in the small sample from H. M. O. Stephenson et al. (2025), as well as the different modelling approaches used to derive the sizes (see Section 4.1 for more discussion).

We compare our rest-frame FUV size measurements to T. Shibuya et al. (2015), focusing on their sample of star-forming galaxies, probed in the UV out to . In the middle panel of Fig. 4, we show their median points for the mass bin, whereas our results are plotted at fixed mass . None the less, our FUV average size kpc is in good agreement with T. Shibuya et al. (2015) at .

Turning to the optical sizes, plotted on the bottom panel of Fig. 4, we find excellent agreement with our measured median size at ( kpc) and the sizes from L. Yang et al. (2025) at this same redshift. Our median measurement lies below those from E. Ward et al. (2024), N. Allen et al. (2025), and T. B. Miller et al. (2025). The overall scatter in the reported median sizes from different observational studies points to a strong effect of the sample selection and modelling done in each work. Furthermore, a single Sérsic modelling approach could be too crude, and is sensitive to the spatial resolution of the data, the sky background, and the contamination from neighbouring objects. The mass ranges probed in the different works are also important. We choose the reference mass as it lies in a well-constrained region of our sample, and fit galaxies with . N. Allen et al. (2025) probe galaxies down to , but as shown on Table 2, our median size is slightly smaller when we fit our full mass range, so it does not help to explain the observed discrepancy. On the other hand, E. Ward et al. (2024) only fits high-mass galaxies , so our reference mass lies at the cut-off of their sample. This is also the case for T. B. Miller et al. (2025), where they instead only fit low-mass galaxies . Finally, L. Yang et al. (2025) use the same mass cut-off as this work and probes a similar range, which can help explain the consistency with our results at .

Across wavelengths, we typically find smaller sizes than simulations. As shown in Fig. 4, our sizes are dex smaller than those reported in W. McClymont et al. (2025b). Notably, they only applied an observability cut based on continuum surface brightness, and so the discrepancy could be caused by observational selections that bias our sample towards more compact objects, missing more extended objects which at fixed SFR will have lower surface brightness. Importantly, the fainter emission at large radii seen in many simulations is typically missed due to the noise and PSF in observations, causing measured sizes to be smaller (P. Punyasheel et al. 2025). This can, at least in part, explain why thesan (X. Shen et al. 2024) and tng50 (L. Costantin et al. 2023) also lie above our median measurements in the FUV and optical. The flares median sizes for galaxies shown on Fig. 4 show good agreement with our measurements in the FUV, but lie well below in the optical. This is mainly due to the turnover in their SMR around , as can be seen on Fig. 3. Overall, their sizes are fairly consistent with our measurements in terms of normalization, but not mass dependence.

On the simulations side, there are various factors which affect the inferred sizes, starting with the chosen method to measure them (e.g. pixel- or aperture-based). Also, based on the strength of the implemented stellar feedback model, galaxies can be more or less bursty, which impacts the distribution of the stars and gas. For example, thesan-zoom is highly bursty. If the burstiness is too strong compared to observations, this could boost the observed sizes due to repeated strong feedback events. The dust content of the simulated galaxies also varies and can have an impact on the measured sizes, which we discuss in detail in Section 4.3.

As shown in Fig. 3, we find that galaxy sizes increase with mass at all probed rest-frame wavelengths (i.e. ), although the slope and normalization vary (Table 2). The sizes have a slightly weaker mass dependence () than the stellar continuum sizes (), although the uncertainties are relatively large. sizes are sensitive to radiative transfer (RT) effects as they probe the extent of the ionization of the interstellar medium (ISM) by young massive stars. The shallower dependence on stellar mass could hence indicate that, at fixed mass, other properties such as SFR, can impact the size.

We now place our measurements in the context of works at both lower and higher redshifts in order to study the redshift evolution of the slope of the SMR. On the top panel of Fig. 5, we plot our slopes at each wavelength, which are constant with redshift, for the and the full sample fits. We represent the uncertainty on the measurements and the redshift span of our sample through the shaded regions.

The shallower slope is consistent with J. Matharu et al. (2024), who find a steeper slope for the rest-frame optical continuum () than for the (). In both cases, we find a slightly steeper slope (consistent within ), which could be attributed to stacking and LSF effects as discussed above. Our measured slopes are also consistent with those reported in H. M. O. Stephenson et al. (2025) at , namely for the optical and for the . Interestingly, W. McClymont et al. (2025b) also find a steeper slope for the continuum compared to . These measurements, which are also consistent with ours within at all comparable wavelengths, suggest that there is an intrinsic, physically driven flattening of the mass dependence. We note that if is tracing non-circular motions such as radial inflows and outflows, this would also affect and perhaps flatten the SMR.

The slopes for the continuum suffer from larger uncertainties, but are consistent with within . We compare our results in the FUV with T. Shibuya et al. (2015). We note that in their work, the growth of sizes is based on UV luminosity and not mass, so this is the slope reported on Fig. 5. This should offer a fair comparison given the linear relation between and at fixed redshift (e.g. E. Daddi et al. 2007). Our measured slope of at lies between the measurements from T. Shibuya et al. (2015) above and below , although consistent within the uncertainties.

In the rest-frame NUV and optical sizes, we find good agreement with other surveys at similar redshifts (T. B. Miller et al. 2025; E. Ward et al. 2024; N. Allen et al. 2025; L. Yang et al. 2025). R. G. Varadaraj et al. (2024) find much flatter slopes () at , which could be due to their ground-based sample selection, which is biased towards luminous systems. In the context of lower redshift measurements (J. Matharu et al. 2024; E. Ward et al. 2024), we find evidence for a constant slope of the SMR. This is consistent with predictions from A. der Wel et al. (2014a) who find a constant slope at . We do not find evidence for a change in slope for the FUV and NUV sizes in our sample, which are in agreement with the optical slope within . This is consistent with the results reported in N. Allen et al. (2025); L. Yang et al. (2025). The lack of wavelength dependence of the slope implies that dust attenuation is weak at these redshifts, or at least has no strong mass dependence. When comparing with the tng50 (L. Costantin et al. 2023), thesan (X. Shen et al. 2024), and flares simulations (W. J. Roper et al. 2022), our observed positive SMR is in tension with the negative or flat slopes reported in these works, as can be seen on Fig. 3. We plot the ‘observed’ SMRs from the simulations, meaning that they are sensitive to dust attenuation effects in the same way as observations, and should offer a more direct comparison. We discuss this further in Section 4.3.

We now focus on the redshift dependence of the SMR, parametrized by in equation (4). We focus on the rest-frame optical and FUV only, from the continuum sizes, since they offer the most comparison samples from the literature. We plot our best-fitting redshift dependence on Fig. 4, as a solid red line, for galaxies at . It is important to note that we are not plotting or discussing the evolution of one galaxy, but rather the evolution of the population of galaxies across cosmic time. Each galaxy will presumably only have a stellar mass of once in its mass assembly history, and we are instead interested in studying whether the average size of a galaxy changes with time. Galaxies with at will have evolved into massive () galaxies at lower redshift, and might even have ceased their star formation.

Starting with , we find a negative but weak redshift dependence , which suffers from large uncertainties due to the small redshift range probed in this work. This implies that the average size of a galaxy remains at kpc at . Comparing to studies at (E. J. Nelson et al. 2016b; D. J. Wilman et al. 2020; J. Matharu et al. 2024), we find that the average size decreases by dex from to . In order for sizes at to evolve from at to kpc at (E. J. Nelson et al. 2016b; J. Matharu et al. 2024), they would have to follow a redshift evolution of , stronger than the one measured in our sample alone, . This points towards a different redshift dependence at high and low redshift, which we discuss more in Section 6.1. None the less, the trend shown on Fig. 4 implies that a galaxy at will be more compact than a galaxy at .

For the stellar continuum sizes, we find a mild positive slope which is consistent with no evolution. Specifically, we infer for the optical sizes, for the FUV, and for NUV. This supports the same conclusion as discussed for the sizes, namely a flattening of the redshift dependence of the SMR at high redshift. Other observational studies predict a range of different slopes. When comparing L. Yang et al. (2025) and E. Ward et al. (2024), neither slope is able to reconcile our measurement at with the sizes from K. V. Nedkova et al. (2021) at . T. B. Miller et al. (2025) also discuss the need for a redshift-dependent slope, although they argue the slope increases at high redshift. They also test a parametrization as a power law in the Hubble constant, , and are still not able to reconcile measurements at low and high redshift with a constant slope .

Differences in the derived redshift evolutions could be attributed to variations in the mass dependence , as well as samples probing different mass and redshift ranges and hence being more or less sensitive to the redshift evolution at . However, it is also important to note that we fit the SMR within our full sample, accounting for the mass and redshift evolution independently. In other works, such as E. Ward et al. (2024) and N. Allen et al. (2025), the redshift evolution is inferred from median sizes (at fixed mass) in different redshift bins. This subtle difference could influence the median sizes inferred, as this latter method does not take into account redshift evolution within a specific bin, and could hence be influenced by bin size, for example. None the less, we also note that our SMR fit (Fig. 3 suffers from relatively large uncertainties, particularly in the redshift dependence parameter , as described above. Samples spanning larger redshift ranges are needed to better constrain this evolution, but this is beyond the scope of this work.

W. McClymont et al. (2025b) argue for the need of a double power-law fit for the sizes in the thesan-zoom simulations, finding a shallow slope at the redshifts probed here. However, exploring a more complex redshift evolution is beyond the scope of this work, as our sample does not span a large enough parameter space. On Fig. 4, we show in grey the evolution as expected for the virial radius . We fix the normalization to the median sizes at measured in this work for  FUV, and optical. Qualitatively, for galaxies, the size evolution does not appear fully consistent with the expected growth of haloes. This suggests that the relation between the sizes of galaxies and their host haloes varies with redshift, as is expected for the relation between stellar and halo masses, traced by the stellar-to-halo-mass (SMHM) relation (P. Behroozi et al. 2019).

Finally, we discuss our results regarding the intrinsic scatter, , of galaxies around the SMR. This scatter sheds light on whether other physical processes influence the sizes of galaxies at fixed stellar mass. We plot our best-fitting measurements on the bottom panel of Fig. 5.

The scatter for all of the continuum sizes is consistent within the uncertainties, , and is larger than its counterpart for the sizes, . Across all wavelengths, the scatter increases when the full sample is fit. This implies that the scatter increases at lower masses, as can be directly seen on Fig. 3. Although this could be due to a number of reasons, it is likely caused by smaller sizes being more difficult to constrain due to the size and complex shape of the PSF. These measurements therefore suffer larger uncertainties. Also, the Sérsic index is more difficult to constrain for smaller systems, and tends to be higher than for larger, more massive galaxies (see Fig. A4). Because the Sérsic index is degenerate with the effective radius, this affects the best-fitting sizes and introduces more variation.

In the bottom panel of Fig. 5, we also plot measurements from other studies at different redshifts. In the NUV and optical, the values from the literature vary between , but the scatter seems to remain at around independently of redshift, consistent with the scatter measured from our fiducial fits. This scatter is however difficult to constrain as it strongly depends on understanding the measurement uncertainties, which can be derived in different ways in the literature. In our work, they are derived in a fully Bayesian framework and marginalized over all of the free parameters. None the less, the relative agreement shown in Fig. 5 suggests that the measured scatter has a distinct physical nature, which we discuss in Section 6.2. In the FUV, T. Shibuya et al. (2015) find a significantly larger scatter from . We also note that their scatter is defined as the standard deviation of the lognormal distribution of sizes at fixed , so we convert it to a decimal (base 10) logarithm to include it in our comparison. They do not fit for the intrinsic scatter but instead account for the measurement errors in quadrature, which could in part explain the discrepancy with our measurements. Importantly, as mentioned, they measure the scatter in the relation, and not the relation, which could introduce systematics.

In numerous studies at lower redshifts, the ratio of to rest-frame UV or optical sizes has been used to motivate the inside–out growth of galaxies (e.g. E. J. Nelson et al. 2016b; D. J. Wilman et al. 2020; J. Matharu et al. 2024). emission is driven by the ionization of gas by young massive stars, through the emission of LyC photons, and hence traces star formation on short time-scales ( Myr), reflective of these stars’ lifetime. The rest-frame UV traces star formation on longer time-scales ( Myr), and finally, the rest-frame optical is a good tracer of the galaxy’s integrated star formation history (SFH). Studying the ratio of galaxy sizes at these different wavelengths provides insight into how star formation has evolved spatially within the galaxy. At higher redshifts, however, galaxies are younger, more gas-rich, and often dominated by recent bursts of star formation. In such systems, strong outshining from young stellar populations and rapid structural evolution may alter the relative sizes measured at different wavelengths, potentially leading to different interpretations than at lower redshifts ().

In the left panel of Fig. 6, we plot the ratio of to NUV sizes as a function of the galaxy’s offset from the main sequence, . Most of our sample (77 per cent) has larger  sizes compared to the NUV, . In terms of star formation time-scales, this implies that younger stars are forming in the outskirts, beyond where the earlier generation had formed. However, following this logic, we should see the same trend when looking at the ratio of FUV to optical sizes. We plot as a function of on the right panel of Fig. 6. As highlighted by the running medians, we find that our sample is scattered around , providing no evidence for steady inside–out growth. Based on Fig. 4, the shallow increase of average sizes at fixed mass, at least at high redshift, would also support this statement. This implies that galaxies could be growing inside–out on average, but both the optical and UV are dominated by the last burst.

This result suggests that the larger sizes, compared to the UV, are probably explained by another phenomenon. The left panel of Fig. 6 shows an increase of above the main sequence, . In fact, we find a mild but significant correlation between these two quantities (, ), and see a clear increase in the running medians from to . We also find evidence for a slight increase below the main sequence, although our sample size drops significantly in that region. When comparing to the thesan-zoom simulations (W. McClymont et al. 2025b), we find good agreement with this observed trend. As discussed in W. McClymont et al. (2025b), the boost in sizes above the main sequence can instead be explained by large amounts of LyC photons produced during bursts of star formation, which increases the radii of their Strömgren sphere (B. Strömgren 1939). The Strömgren sphere describes the ionizing front of the LyC photons, which becomes the radial extent of the induced emission. Hence, the sizes do not directly trace their source of ionization, i.e. the young massive stars, but rather reflect the intensity of the ongoing star formation. During a burst, where is high, the sizes will then naturally be at their highest. Given the lack of boost of NUV sizes compared to the optical, this phenomenon is a more likely explanation for the larger sizes than inside–out growth. This is supported by observations finding that star formation is bursty at high redshift (e.g. T. J. Looser et al. 2025; S. Tacchella et al. 2023; C. Simmonds et al. 2025; W. M. Baker et al. 2025b), which promotes galaxy growth through central starbursts (K. El-Badry et al. 2016; S. Tacchella et al. 2016b; W. McClymont et al. 2025b) rather than a star-forming disc. We find evidence for this in our work, with the observed increase of with is in part driven by the decrease of UV sizes during bursts (Fig. 7). The results are also consistent with medium-band imaging observations, which show that the relative size of emission line regions compared to the continuum scales with the ionizing photon production rate (Y. Zhu et al. 2025).

To qualitatively explore the validity of our Strömgren sphere interpretation, we estimate the Strömgren radii for the galaxies in our sample using a simple approximation. The Strömgren radius is empirically defined as , where is the ionizing flux, n is the hydrogen number density, and is the total recombination rate. Assuming a Case B recombination rate for an ISM temperature of K, and an average density of 10 cm, we can estimate using the integrated flux, scaled by the redshift to obtain the intrinsic luminosity. We find a median kpc, which is comparable to our measured radii, meaning our interpretation is plausible.

Our finding of large to NUV size ratios at is consistent with measurements from Villanueva et al. (in preparation) at , where they model and subtract the rest-optical continuum from the F356W photometry, by using non-contaminated neighbouring medium band filters, to obtain maps. They attribute the large sizes to a combination of vigorous star formation and leaking ionizing photons. Interestingly, H. M. O. Stephenson et al. (2025) do not find that the nebular () sizes are larger than the stellar continuum (UV or optical) sizes on average, as they are instead found to be typically smaller than the stellar continuum traced by the rest-frame R-band (F444W), with a median . The small sample size and large scatter in the measurements, which could be in part caused by the low stellar masses probed (), may be responsible for the discrepancy with our results. Furthermore, their selection is based on narrow-band H-alpha emission, which could cause a bias towards galaxies with centrally concentrated profiles. Discrepancies could also be caused by the different modelling approaches, where H. M. O. Stephenson et al. (2025) fit the narrow-band imaging, containing the emission, with galfit, and we forward-model the emission line in the grism data.

A key factor to consider when studying the dependence of multiwavelength size ratios with (Fig. 6) is the behaviour of the individual sizes themselves. We investigate the rest-frame FUV sizes, as they are expected to probe star formation on the shortest time-scales. As shown in Fig. 7, we find a clear decrease of with , suggesting that galaxies go through a compaction as they move above the main sequence. However, selection effects lead to a form of degeneracy between stellar mass and offset from the main sequence in our sample (Fig. 1), which means the distinct effects of these two quantities become intertwined. We attempt to break this degeneracy using our Bayesian fitting approach with emcee.

As described in Section 3, we fit our SMR by also accounting for a redshift dependence (equation 4). To quantify the dependence of sizes on both stellar mass and , we compute a new fit of the SMR where we further include a power-law dependence on by adding a term in equation (4), since is already in log-space (equation 5). We show the posterior distributions for , , (equation 4) and the new parameter on Fig. A3. Compared to the fit without the dependence (Table 2), we find that the mass dependence is slightly weaker in this new fit, from to , although consistent within the uncertainties. More importantly, we find a weak positive, non-zero dependence of on , . This implies that within our data, the FUV sizes decrease above the main sequence independently of their stellar mass or redshift.

We plot the best-fitting relation at fixed redshift () and for a range of stellar masses () on Fig. 7. As expected, we find the intrinsic dependence to be weaker than the trend traced by the medians, meaning some of the observed decrease is driven by stellar mass. Larger, more complete samples are needed to better quantify the dependence of sizes on and . Our running medians provide a form of upper limit on the strength of the relation, since they implicitly assume that is the only factor. When we account for mass in our sample, we find that the dependence on is weaker, but this could be underestimated due to a biased sample at low masses.

The compaction of galaxies as they move up the main sequence has been studied in cosmological simulations (A. Zolotov et al. 2015; S. Tacchella et al. 2016b; S. Lapiner et al. 2023; W. McClymont et al. 2025b). This state is characterised by compact, gas-rich cores of star formation, which occur when the galaxy reaches its peak star formation. This strong, centrally concentrated star formation then leads to temporary inside–out quenching where feedback drives the gas from the core (S. Tacchella et al. 2016b; W. McClymont et al. 2025a,2025b). These oscillations between compaction and quenching predict that should influence the observed sizes of galaxies independent of their mass, with galaxies above the main sequence being seen preferentially in compaction phases. In this work, we find evidence of this trend, despite significant scatter and uncertainties. Studying where quiescent galaxies, which lie well below the main sequence, lie on these trends would be highly informative but is beyond the scope of this work since we are limited by the need for detections. In this work, we only probe systems that are in the process of being temporarily quenched, meaning they have falling SFHs.

Even though most observations find positive SMRs across redshifts (Section 3.3), cosmological hydrodynamical simulations have struggled to reproduce this result as they find flat or negative slopes at high redshifts (M. A. Marshall et al. 2022; W. J. Roper et al. 2022; L. Costantin et al. 2023; X. Shen et al. 2024), as shown in Fig. 3. This discrepancy is often attributed to the effects of dust attenuation, which increases with stellar mass and is stronger in the central regions of the galaxy (X. Wu et al. 2020; M. A. Marshall et al. 2022; G. Popping & C. Péroux 2022; W. J. Roper et al. 2022; L. Costantin et al. 2023). When they take into account dust, they are able to increase their slope because the dust boosts the sizes of massive galaxies by obscuring the central regions. P. Punyasheel et al. (2025) show that when accounting for observation effects, such as noise and PSF, using the FLARES simulations, the slope then even becomes positive.

Recently, W. McClymont et al. (2025b) found an intrinsically positive SMR in the thesan-zoom simulations, and argued that this was recovered due to the bursty nature of star formation in those simulations allowing for regular compaction and expansion as galaxies build-up mass. This aligns with the findings of W. J. Roper et al. (2023), who claim that the negative SMR at high redshift in large-volume cosmological simulations with an effective equation-of-state ISM model is due to runaway gas collapse down to the gravitational softening length in the most massive, highly star-forming galaxies at high redshift in these models. This raises the interesting prospect that we can constrain models of galaxy formation by observationally measuring the impact of dust on high-redshift galaxy sizes.

In order to investigate whether our recovered slopes of could be due, at least in part, to dust attenuation, we show the ratio of FUV to NUV sizes in Fig. 8 as a function of stellar mass. We colour-code our galaxies by the measured dust attenuation, extracted from our SED modelling results, which is represented by the optical depth . These two rest-frame wavelength ranges probe similar physical processes at high redshift, but the bluer wavelengths (FUV) are more affected by dust attenuation. If our galaxies were affected by central dust obscuration, we would expect to see a positive value of , as the FUV radii would appear larger. This effect should be more extreme for higher mass galaxies, which typically have more dust. However, as highlighted by the running medians on Fig. 8, our sample is scattered around , showing no systematic enhancement of FUV sizes. The scatter around is consistent with measurement errors and the intrinsic scatter measured around the SMR (Fig. 5, bottom panel). The measured attenuations from SED modelling with prospector show that most of our galaxies have relatively low dust attenuation (), and, as expected, the higher mass galaxies have typically higher values of . However, Fig. 8 suggests that the dust attenuation in our galaxies is homogeneous, and not centrally concentrated, and therefore does not introduce biases in our measured sizes.

We compare our measurements to the FUV-to-NUV ratio measured in the flares simulations (W. J. Roper et al. 2022). As can be seen on Fig. 8, there is a trend of increasing with stellar mass, although our sample is scarce at the high masses where this becomes significant in flares. None the less, this boosting at high masses is consistent with the galaxies having significant central dust obscuration. It highlights that the intrinsic SMRs would be strongly affected by the increasing dust attenuation for higher mass galaxies. We qualitatively also compare to the thesan simulations (X. Shen et al. 2024) by plotting the ratio of observed UV to intrinsic UV (no dust) sizes. This is not an apples-to-apples comparison, but it illustrates the effects of dust attenuation well. The lack of such a trend in our data suggests that our measured SMR is intrinsically positive and not affected by dust. These comparable NUV and FUV sizes are also found in Villanueva et al. (in preparation) at similar redshifts ().

We now investigate how the evolution of sizes at different wavelengths is related to the kinematic evolution of galaxies, namely their ionized gas rotational support and intrinsic velocity dispersion . As described in Section 2.6, these kinematic quantities are obtained through kinematic modelling with geko and presented in A. L. Danhaive et al. (2025).

In Fig. 9, we study how the ratio evolves with the rotational support of the ionized gas, . Specifically, we want to relate the rotational support to the movement of galaxies around the main sequence. traces whether galaxies are above or below the star-forming main sequence, on a 10 Myr time-scale, but does not contain information on their direction, i.e. . Although this is a difficult quantity to constrain, we can approximate it by measuring the ratio of the SFR averaged over the last 10 Myr and over the last 30 Myr (). If a galaxy is moving above the main sequence, then its SFH should be rising, hence . Conversely, if it is moving below, then we should observe . We extract these two quantities from our SED modelling (Section 2.4) and colour-code galaxies on Fig. 9 based on whether they have rising (purple) or falling (yellow) SFHs. We note that this cut splits our samples roughly in half, with 55 per cent of galaxies having rising SFHs, which is expected since galaxies should not be found preferentially in either state.

Interestingly, we find that only correlates with for galaxies with rising SFHs (, ), with no trend measured for galaxies with falling SFHs (, ). This can also be seen through the running medians, which are rising and constant, respectively. If we adopt the picture painted in Fig. 6, where galaxies do not grow inside–out but instead through central starbursts, then Fig. 9 can be interpreted as accreting gas being illuminated by the LyC ionizing photons emitted from the central regions. As the sizes are boosted, they reach further out and could be probing the gas that is fuelling the ongoing star formation. The increase of would then suggest this gas is not infalling radially but rather following circular motions. When galaxies start moving back down towards and below the MS, it implies that star formation is being shut off in the centre, most likely due to outflows and feedback (S. Tacchella et al. 2016a; W. McClymont et al. 2025a,2025b). On longer time-scales, these processes introduce turbulence in the gas (Y. Mai et al. 2024; A. L. Danhaive et al. 2025) and drive the decrease of . Also, with the decrease of the size of the region, we become more sensitive to the central regions of the galaxy where gas is not following circular motions. We discuss this further in Section 6.

A key question when measuring gas kinematics, and even more so ionized gas kinematics, is about what exactly we are probing and measuring. For example, mergers can often imitate the kinematic signatures of rotating discs (B. E. Robertson & J. S. Bullock 2008; R. C. Simons et al. 2019), and the rotational support can vary by orders of magnitude depending on the gas tracer used for the measurement (M. Kohandel et al. 2020; F. Rizzo et al. 2024). At lower spatial resolutions, the effects of outflows and inflows can also impact kinematic measurements, making it difficult to distinguish the origins of the velocities measured.

In Section 5.1, we argue that the trends observed in Fig. 9 could imply that galaxies build up rotational support, as they accrete gas to fuel star formation. As the sizes decrease after the starburst, the rotational support decreases due to the smaller scale probed and the turbulence introduced by processes such as feedback from star formation. In this view, the rotational support at high redshift, as measured through , is more of a probe of the baryon cycle than the build-up of discs. Because of the observational bias promoting high SFRs at low masses in our sample (Fig. 1), such an interpretation could explain the unexpected decreasing trend observed between and in A. L. Danhaive et al. (2025).

Many of the statistical claims of discs in the early universe are based on morphological studies (L. Ferreira et al. 2022; B. E. Robertson et al. 2023), which often use projected axis ratios and Sérsic indices to infer the disc-like nature of galaxies in the absence of kinematic measurements. As such, galaxies with small values of and/or Sérsic indices close to one are often categorized as discs. E. J. Nelson et al. (2023) and J. L. Gibson et al. (2024) discuss the potential disc-like nature of red ultra-flattened objects uncovered by JWST, and further raise the question of how well the morphology traces the underlying kinematics.

To investigate this further, we plot our measured UV axis ratios as a function of the rotational support in Fig. 10. Across values of , we find a diversity of Sérsic indices n. If the elongated systems in our sample, with , were predominantly discs, we would expect them to have higher values of than those with , which are more spheroidal in nature. However, we instead see the large scatter of values of at fixed . This scatter is in part expected as some objects with rotational support could have high values of simply due to being observed face-on. The running medians and the percentiles shown on Fig. 10 are consistent with a flat trend. We also do not find any trends between rotational support and Sérsic indices measured both in the UV and the . This could in part be caused by biases and uncertainties in the measurement of Sérsic indices at high redshift, where galaxies are smaller both intrinsically and on the sky. Physically, the lack of a visible correlation could be driven by the turbulent and bursty nature of galaxy evolution at high redshift, causing the morphology and kinematics of galaxies to change on short time-scales and not settle along a stable relation.

We find that not all elongated systems are disc-like in nature, but instead that 44 per cent of  systems are flattened systems, which are not rotationally supported. We also plot the distribution of b/a values on Fig. 10, and find that it is not uniform but instead peaking around . Whereas for a population of discs we would expect a flat distribution, such a peaked shape implies the contribution of spheroids and morphologically prolate systems (see e.g. B. E. Robertson et al. 2023; V. Pandya et al. 2024). This result is in agreement with studies such as J. Vega-Ferrero et al. (2024), who use the tng50 simulations to show that visual and/or morphological classification of galaxies overestimates the fraction of disc-like systems, mislabelling systems which actually have low angular moment and non-oblate structures. V. Pandya et al. (2024) also discuss the difficulty in constraining 3D shapes from 2D projections, which can lead to misclassifications. More importantly, these elongated non-rotating systems could be the result of filamentary accretion along the cosmic web, with the baryons tracing prolate DM distributions (D. Ceverino, J. Primack & A. Dekel 2015). We also note that although these systems are not rotationally supported, some of them have resolved rotation which is simply secondary with respect to the pressure support. These systems could also be rotating around their minor axis (‘kinematically’ prolate), although this has never been observed.

The lack of low systems with high is most likely due to the observational limitation of only measuring projected velocities, due to galaxy inclination. The true rotational velocity is obtained by de-projecting the observed velocity, , so low values of i (i.e. high values of ) would boost . This results in large uncertainties, as can be seen in Fig. 10, but can explain the observed bias.

Analytical theoretical models of galaxy formation provide predictions of the SMR, along with its redshift dependence and intrinsic scatter. In Section 3.3, we show that our measurements of the slope of the SMR are consistent with those found down to . As discussed in S. Shen et al. (2003), the slope of SMR can be linked to the ratio of stellar to halo mass, and how it depends on the halo mass. Assuming a simple disc model (S. M. Fall & G. Efstathiou 1980; H. J. Mo et al. 1998), a slope of implies that this ratio is not constant, but rather follows a power law. In fact, assuming and , then where . The lack of redshift evolution in the slope directly implies that this stellar-to-halo mass relation holds out to redshifts of at least , beyond which sample sizes decrease, leading to large uncertainties (see Fig. 5).

However, it is important to consider the stellar mass range probed. L. Mowla et al. (2019) find that the slope of the SMR is actually mass dependent. In fact, when fitting the SMR with a double power-law, they find that the pivot mass evolves from to from to . Specifically, the evolution of these pivot masses seems to trace the pivot masses of the SMHM relation from P. Behroozi et al. (2019). In our work, we do not probe high enough masses to probe this double power-law behaviour at , but it is an interesting scope for future study.

In terms of the redshift dependence, we show in Section 3.4 that the measured flat redshift dependence found in this work, as well as the evolution of the median size at between and , suggest that changes with cosmic time. Specifically, we find evidence for a flattening in the growth of sizes at high redshift, although larger samples, spanning larger redshift ranges, are needed to explore this further. Such a flattening could be attributed to the lack of stable gaseous discs at high redshift (A. L. Danhaive et al. 2025), driven by bursty star formation, mergers, and high gas fractions, amongst others. In fact, as galaxies are able to settle into thick discs at cosmic noon and then thin discs at , their sizes at fixed mass would increase as the radial extent of the disc increases. On the other hand, at high redshift, galaxies are not able to grow smoothly inside–out due to central starbursts. Overall, a change of with redshift would imply that different underlying mechanisms influence the sizes of galaxies at different redshifts. None the less, studies of sizes across cosmic time are consistent with the growth of galaxy sizes being driven primarily by the growth of the virial radius of the underlying halo.

Turning to the intrinsic scatter around the SMR, cosmological simulations and theoretical models provide a good framework for understanding its origin. In fact, this scatter is often linked to the distribution of spin parameters of DM haloes at fixed mass (H. J. Mo et al. 1998; A. der Wel et al. 2014a). It has been measured within simulations to be around dex (e.g. A. Burkert et al. 2016), which is consistent with our best-fitting values. We find that the offset of galaxies from the SMR follows a lognormal distribution (see Fig. A5), which is also consistent with this interpretation. This suggests that, at least at the masses probed by our fiducial fits (), the scatter around the SMR can be attributed to the variation in halo spins, although detailed analysis is required to test this hypothesis (see also M. Danovich et al. 2012, 2015).

Additionally, the bursty nature of star formation at high redshift could increase the scatter of sizes at fixed mass, particularly at wavelengths more sensitive to recent star formation (i.e. FUV and NUV). Central bursts of star formation are not necessarily spatially correlated to the size or mass of the galaxy, and would have a bigger impact for smaller galaxies, which could help explain the increase in scatter when including galaxies with . The region of ionized gas is predicted to extend beyond the direct star-forming regions (see Section 4.1), and instead is more dependent on the galaxy’s SFR which describes the amount of ionizing photons produced. C. Simmonds et al. (2025) conduct a detailed study of the SFMS and find that the scatter related to short-term burstiness, the difference between the intrinsic scatter of the SFMS between a look-back time of 10 and 50 Myr, is of the order of dex in our mass and redshift range. This scatter is consistent with the intrinsic scatter around the SMR, suggesting a potential link between the distribution of sizes at fixed mass and the burstiness of star formation.

In Sections 3, 4.1, and 5, we study the relation between sizes at different rest-frame wavelengths and tracers of ongoing star formation and ionized gas kinematics. Specifically, we have focused on how these different properties vary as galaxies move along the MS ridge-line, and what that can teach us about galaxy growth at high redshift. We discuss here a coherent interpretation of our results, although larger sample sizes and joint studies with simulations are needed to further test its validity. In the local Universe, and up to cosmic noon, weak inside–out growth has been observed through the increase of equivalent width or sSFR with radius (S. Tacchella et al. 2015; E. J. Nelson et al. 2016b; F. Belfiore et al. 2018; L. Morselli et al. 2019), implying that new stars are forming preferentially in extended star-forming discs. This stable growth is typically observed in galaxies with . At high redshift , galaxies have been shown to have disturbed kinematics (A. L. Danhaive et al. 2025), bursty SFHs (e.g. T. J. Looser et al. 2024; J. Trussler et al. 2025; C. Witten et al. 2025; W. M. Baker et al. 2025b), and high gas fractions (L. J. Tacconi et al. 2018; L. J. Tacconi, R. Genzel & A. Sternberg 2020; E. Parlanti et al. 2023; A. Graaff et al. 2024a; A. Graaff, A. Pillepich & H.-W. Rix 2024b). As such, we do not expect them to grow in the same way as more massive galaxies at lower redshifts. It is important to note that some galaxies at high redshift have thick discs (; A. L. Danhaive et al. 2025), but their low values of raise the question of whether these discs are long-lived, or whether they are periodically destroyed.

As shown in the right panel of Fig. 6, we find that the UV traces the optical at , similar to what has been found in other works at high redshift (Y. Ono et al. 2024; K. Ormerod et al. 2024; N. Allen et al. 2025; L. Yang et al. 2025). We also show in Fig. 8 that this is unlikely to be caused by central dust obscuration of the UV continuum. Together, these results do not favour inside–out growth. Furthermore, inside–out growth is understood as the formation of a central bulge, followed by the development of an extended star-forming disc. However, as discussed in Section 5.2, it is unlikely that star-forming galaxies are forming long-lived discs at in the mass range probed in this work ().

Instead, our work suggests an evolution driven by the movement of galaxies about the SFMS ridge-line. In Fig. 11, we show a schematic view of the evolution of and UV sizes as galaxies go through a burst of star formation. As galaxies accrete gas towards the central regions, it forms stars. Their LyC photon emission ionizes the gas, forming an halo which expands as the star formation increases. Therefore, as shown on Fig. 6, the higher a galaxy is above the main sequence, the larger its -to-NUV size ratio is. The increase of measured as increases above the main sequence suggests that the illuminated gas is following circular motions as it is accreting. Once stellar feedback starts to kick in, it injects turbulence into the gas and drives it from the central regions through outflows. These effects could explain the lack of trend as galaxies move down the main sequence (Fig. 9), as the gas undergoes less ordered motions. Furthermore, depending on the mass of the galaxies and the amount of star formation, the size of the region could be too small to probe the inflow of gas on large scales, which instead could maintain rotational support. Within this framework, acts as a tracer of the baryon cycle within galaxies taking place in the central regions. Cold gas (i.e. CO or [CII]), on the other hand, typically probes larger scales which are less affected by short-term variations in star formation.

Another plausible scenario is that outflows are driving the extended sizes. In fact, studies of [CII] emission with the Atacama Large Millimeter Array (ALMA) also find larger [CII] sizes compared to UV, consistent with gas on large scales being illuminated by young stars at the centres of galaxies (S. Carniani et al. 2018; S. Fujimoto et al. 2020). These [CII] detections could be tracing similar gas to the . However, this large-scale gas is typically associated to outflowing, and not inflowing, gas, based on the elevated velocity dispersions. Identifying and characterizing outflows at high-redshift is challenging, particularly if they have low velocities, but we expect strong central star formation to result in outflows (S. Carniani et al. 2024) which expel the gas from the central regions (W. McClymont et al. 2025a,2025b). In this work, we favour the inflow scenario mainly based on the fact that the trend of increasing size with is only seen for galaxies with rising SFHs. In the outflow-driven case (Fig. 9), we would expect the outflows to predominantly kick in when the SFR starts to decline. Observing a trend for galaxies with falling SFHs could have therefore supported the outflow scenario. None the less, outflows could still be contributing to the observed trends and high-resolution IFU data for large samples are needed to distinguish these two scenarios.

There are some caveats of this work that are important to mention in the context of this discussion. As shown in Fig. 1, we suffer from a bias towards highly star-forming systems at low masses. This is the case for most observational studies, and implies a correlation between and . We run a partial correlation coefficient (PCC) analysis for the ratio and find that neither quantity ( or ) can fully account for the other’s influence on this ratio. However, we lack the appropriate sample sizes needed for PCC analysis, as can be seen by the bias in the random variable. Aside from the PCC analysis, where we expect sizes to trace stellar mass through the SMR, we do not expect the ratio of sizes to be as sensitive to stellar mass, but rather to physical processes like star formation.

As discussed in Section 3.2, single Sérsic fits could be too simplistic, and hence introduce some difficult-to-quantify biases. This is, for example, the case for galaxies with a bulge and disc, which typically have different Sérsic indices for each component. Although should capture sizes to first order, studies of clumpiness and the use of multicomponent modelling are needed to test the validity of single Sérsic fits.

There are also uncertainties relating to the kinematic measurements, which could impact our conclusions. They are discussed in detail in A. L. Danhaive et al. (2025), but the main ones are caused by resolution effects and the lack of full 3D data. It is difficult, given grism data, to constrain non-circular motions such as inflows and outflows of gas, which we know will be playing a role in galaxies at high redshift (S. Carniani et al. 2024). Quantifying their potential imprint on grism data requires careful forward modelling from cosmological simulations, which is beyond the scope of this work.

Finally, as highlighted in Fig. 1, our sample’s selection suffers from strong degeneracies between stellar mass and offset from the main sequence. Although we attempt to separate the effects of these two quantities in this work, larger samples would be useful to statistically quantify their relative importance further.