# Stellar feedback in the star formation-gas density relation: Comparison between simulations and observations

P. Suin, A. Zavagno, T. Colman, P. Hennebelle, A. Verliat, D. Russeil

► To cite this version:

P. Suin, A. Zavagno, T. Colman, P. Hennebelle, A. Verliat, et al. Stellar feedback in the star formation-gas density relation: Comparison between simulations and observations. Astronomy & Astrophysics, 2024, 682, 10.1051/0004-6361/202347527. insu-04479007

HAL Id: insu-04479007

https://insu.hal.science/insu-04479007

Submitted on 8 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stellar feedback in the star formation–gas density relation: Comparison between simulations and observations

P. Suin<fnref n="1" />, A. Zavagno<fnref n="1" />,<fnref n="2" />, T. Colman<fnref n="3" />, P. Hennebelle<fnref n="3" />, A. Verliat<fnref n="3" />, and D. Russeil<fnref n="1" />

Received 21 July 2023 / Accepted 24 November 2023

## ABSTRACT

**Context.** The impact of stellar feedback on the Kennicutt–Schmidt (KS) law, which relates the star formation rate (SFR) to the surface gas density, is a topic of ongoing debate. The interpretation of high-resolution observations of individual clouds is challenging due to the various processes at play simultaneously and inherent biases. Therefore, a numerical investigation is necessary to understand the role of stellar feedback and identify observable signatures.

**Aims.** In this study we investigate the impact of stellar feedback on the KS law, aiming to identify distinct signatures that can be observed and analysed. By employing magnetohydrodynamic simulations of an isolated cloud, we specifically isolate the effects of high-mass star radiation feedback and protostellar jets. High-resolution numerical simulations are a valuable tool for isolating the impact of stellar feedback on the star formation process, while also allowing us to assess how observational biases may affect the derived relation.

**Methods.** We used high-resolution ($<$0.01 pc) magnetohydrodynamic numerical simulations of a $10^4 M_\odot$ cloud and followed its evolution under different feedback prescriptions. The set of simulations contained four types of feedback: one with only protostellar jets, one with ionising radiation from massive stars ($>$8 $M_\odot$), one with the combination of the two, and one without any stellar feedback. In order to compare these simulations with the existing observational results, we analysed their evolution by adopting the same techniques applied in the observational studies. Then, we simulated how the same analyses would change if the data were affected by typical observational biases: counting young stellar objects (YSO) to estimate the SFR, the limited resolution for the column density maps, and a sensitivity threshold for detecting faint embedded YSOs.

**Results.** Our analysis reveals that the presence of stellar feedback strongly influences the shape of the KS relation and the star formation efficiency per free-fall time ($\epsilon_{ff}$). The impact of feedback on the relation is primarily governed by its influence on the cloud's structure. Additionally, the evolution of $\epsilon_{ff}$ throughout the star formation event suggests that variations in this quantity can mask the impact of feedback in observational studies that do not account for the evolutionary stage of the clouds. Although the $\epsilon_{ff}$ measured in our clouds is higher than what is usually observed in real clouds, upon applying prescriptions to mimic observational biases we recover a good agreement with the expected values. From that, we can infer that observations tend to underestimate the total SFR. Moreover, this likely indicates that the physics included in our simulations is sufficient to reproduce the basic mechanisms that contribute to setting $\epsilon_{ff}$.

**Conclusions.** We demonstrate the interest of employing numerical simulations to address the impact of early feedback on star formation laws and to correctly interpret observational data. This study will be extended to other types of molecular clouds and ionising stars, sampling different feedback strengths, to fully characterise the impact of H II regions on star formation.

**Key words.** methods: numerical – stars: formation – ISM: clouds – HII regions – ISM: structure

## 1. Introduction

The existence of a power-law relation between the star formation rate (SFR) per unit area and the gas surface density has long been established through observations of external galaxies (Schmidt 1959, 1963; Kennicutt 1989). The relation, known as the Kennicutt–Schmidt (KS) law, remarkably holds on scales spanning from entire galaxies to individual molecular clouds.

The dichotomy observed between starburst and non-starburst galaxies (Gao & Solomon 2004; Bigiel et al. 2008; Krumholz et al. 2012) might be a hint that stellar feedback plays a role in setting the star formation law (Ostriker & Shetty 2011), and Kennicutt & De Los Reyes (2021) suggested that this bimodal (or multi-modal) relation could originate from changes in the

small-scale structure of the molecular interstellar medium. Indeed, even nearby galaxies show variations on 1.5 kpc scales when observed at higher angular resolution (Sun et al. 2023).

In effect, high-mass stars' feedback, such as photo-ionising radiation, stellar winds, and supernova explosions, greatly impacts the surrounding environment (Chevance et al. 2022), but how it affects the star formation process is still not completely understood. For example, H II regions are capable of blowing out the molecular cloud from which they were born. Yet, their expansion locally compresses the gas, enhancing its density and directly impacting the properties of newly forming clumps (Zhang et al. 2020, 2021). These events can affect scales from tens of parsecs down to sub-parsec scales, possibly altering the local star formation properties (Palmeirim et al. 2017; Mazumdar et al. 2021), or even promoting the formation of

<footnote n="1">¹ Aix Marseille Univ, CNRS, CNES, LAM Marseille, France
e-mail: paolo.suin@lam.fr</footnote>

<footnote n="2">² Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>

<footnote n="3">³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris-Diderot Sorbonne Paris-Cité, 91191 Gif-sur-Yvette, France</footnote>
other massive stars by inhibiting gas fragmentation (Hennebelle et al. 2022).

Recent observations from the James Webb Space Telescope of the NGC 628 galaxy (Barnes et al. 2023; Mayya et al. 2023) have shed light on the importance of understanding the role of early and smaller-scale high-mass stellar feedback for a reliable description of star formation laws. However, effects acting on such small scales are invisible when averaged on galactic scales but may be significant when considering singular molecular clouds. For this reason, observational studies have looked at the KS relation in several low- and high-mass molecular clouds (Gutermuth et al. 2011; Willis et al. 2015; Pokhrel et al. 2021; Bieging & Kong 2022), reporting a steepening of the exponent $\beta$ when passing from the largest galactic scales ($\beta \approx 1.5$; Kennicutt & Evans 2012; Kennicutt & De Los Reyes 2021) to the smallest ones ($\beta \approx 2$; Evans et al. 2014).

Many theories attempt to explain the relation in all its shapes, but verifying or even constraining them represents a challenge for observational studies. Aside from resolution and sensitivity biases, the use of different tracers to derive the local physical conditions, and the multitude of processes ongoing simultaneously, make the comparison with the idealised cases extremely complicated. Moreover, these studies often rely on assumptions to compensate for the lack of 3D information and uncertainties in theoretical models, which hinders high-precision measurements of the real SFR-gas density relation.

For instance, one of the most popular theories is the model proposed by Krumholz & McKee (2005). Under the assumption that the main force governing the process is gravity, the model proposes a linear dependence between $\dot{\Sigma}_*$ and $\Sigma_{gas}$, with the typical timescale given by the free-fall time, $t_{ff}$, of the region. Mathematically, this can be written as

$$ \dot{\Sigma}_* = \epsilon_{ff} \frac{\Sigma_{gas}}{t_{ff}}, \qquad (1) $$

where $\epsilon_{ff}$, the star formation efficiency per free-fall time, is the only parameter required. All mechanisms other than gravity, feedback included, act to modify this coefficient. Despite its simplicity, the model yielded robust results (Krumholz et al. 2012; Usero et al. 2015; Utomo et al. 2018), opening the path to more refined theories that account for diverse local environmental properties. In particular, models incorporating the concept of multi-free-fall time are particularly promising (Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart 2018). Unlike the original model, which relies on average density, these newer theories infer a density distribution to more accurately describe the gaseous environment. These models not only align more closely with observational data, but they also help explain the higher rates of star formation observed in starburst galaxies (Federrath 2013). Nevertheless, although elegant from an analytical perspective, using the local free-fall time represents a major challenge for astronomers. Specifically, obtaining this quantity typically requires resorting to strong assumptions to estimate the average volume density, such as the uniform-sphere approximation.

Even deriving the local SFR is complicated. In recent works conducted on nearby Galactic clouds, this is inferred by counting the number of detected young stellar objects (YSOs), but their mass and age have to be assumed (e.g. Lombardi et al. 2013; Lada et al. 2017). Additionally, in the densest regions of clouds, the count of YSOs is bound by the instruments' sensitivity (Megeath et al. 2016). All these problems might explain why

some studies claim that this parameter remains constant at all column densities (Pokhrel et al. 2021), while others suggest that the model cannot provide a satisfactory description of the process in single molecular clouds (Evans et al. 2014; Lee et al. 2016; Bieging & Kong 2022).

Part of the problem could be that the molecular clouds used in these studies vary in mass, structure, star formation activity, and evolutionary stage. Some contain high-mass stars, while others are known to be low-mass star-forming regions. Encompassing all these effects simultaneously, the individual impacts on the star formation laws blur together. If stellar feedback actually plays a role, accurately separating the clouds according to their physical condition could solve the tension.

To overcome these challenges, we employed high-resolution numerical simulations to study the possible impact of different types of stellar feedback on the KS relation across the temporal evolution of one cloud. Using simulations allowed us to investigate the time evolution without having to deal with many of the difficulties associated with observations. Our simulation suite consisted of two runs where either protostellar jets or ionising radiation is modelled, one run in which both feedback mechanisms are included, and one simulation without any stellar feedback. In this way, we could analyse the effect of each of these feedback mechanisms separately.

The primary objectives of our work are twofold. First, we aimed to characterise the impact of stellar feedback on the star formation laws at the cloud scale and to disentangle the temporal evolution of these effects. Second, we sought to identify observational signatures associated with these phenomena. To achieve this, we analysed the simulations using observational techniques and compared our findings with recent observational results. Additionally, we considered the biases that may affect these procedures, focusing on the impact of YSO counting on the estimate of the SFR, the low resolution of column density maps, and the sensitivity limit for detecting faint and extinct stars. The purpose of our investigation is to improve our understanding of the intricate relationship between feedback, the evolutionary stage, and the observed properties of star-forming clouds.

The paper is organised as follows. In Sect. 2, we describe the simulations, and in Sect. 3 we explain the methods adopted to derive the KS relation and $\epsilon_{ff}$. Section 4 contains the results of the analysis, and in Sect. 5 we investigate the possible impact of observational biases on the observed relations. Finally, in Sect. 6 we contextualise the obtained results from a physical point of view. Conclusions are given in Sect. 7.

## 2. Numerical simulations

In our study, we analysed the simulations presented in Verliat et al. (2022). This set of four runs followed the evolution of the same isolated $10^4 M_{\odot}$ cloud under the influence of diverse feedback prescriptions: one with both ionising radiation and protostellar jets, one with only ionising radiation, one with only protostellar jets, and one without any feedback. For simplicity, we refer to these runs as HIIR+PSJ, HIIR, PSJ, and NF. Thanks to the wide range of physical processes implemented and the high resolution achieved (of the order of thousands of AU), these simulations provide an ideal test-bed for investigating the effect of stellar feedback on the star formation properties and for testing the observational techniques employed in the field. Here we briefly describe the global properties and the setup of the model, referring to Verliat et al. (2022) for a complete description of how the various physical processes are implemented.
## 2.1. Code and numerical parameters

The magnetohydrodynamic simulations were carried out using the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). The domain was set up as a cubic box of $L = 30.4$ pc, with open boundaries and an initial resolution of $128^3$ cells. Adopting five AMR levels, the maximum resolution achieved in each simulation is $7.4 \times 10^{-3}$ pc ($1.5 \times 10^3$ AU). The refinement criterion assures that the local Jeans length is always resolved with at least 40 cells. When the density of a cell at the highest AMR level exceeds $10^7$ cm$^{-3}$, a sink particle is created (Bleuler & Teyssier 2014). These sinks interact with the environment by accreting material from their surroundings and ejecting 1/3 of it under the form of protostellar jets. Each time the sinks accrete $120 M_\odot$ of gas overall, an ionising star is created inside one of them. The stellar mass is drawn from a Salpeter (1955) distribution between 8 and $120 M_\odot$. Then, the parametric mass-UV luminosity relation given in Vacca et al. (1996) sets its ionising flux. The sequence of stellar masses extracted had been preserved between simulations HIIR and HIIR+PSJ so that they formed the same stars throughout the cloud’s evolution.

## 2.2. Initial conditions

Initially, the cloud was modelled as an approximated $10^4 M_\odot$ Bonnor–Ebert sphere, with a diameter of 15.2 pc and a scale radius $r_0 = 2.5$ pc. Its density was distributed according to

$$n(r) = \frac{n_0}{1 + \left(\frac{r}{r_0}\right)^2}, \qquad (2)$$

with $r$ the distance from the box centre, and the central density $n_0 = 800$ cm$^{-3}$. The rest of the cube is filled with a uniform density of $8$ cm$^{-3}$ up to $r = L$, and the density in the remaining corner regions is set to $1$ cm$^{-3}$.

The gas thermal behaviour was parameterised with a cooling function that takes into account several chemical processes (see Audit & Hennebelle 2005 for a detailed description). This set the initial temperature of the medium at 10 K in the densest parts. The magnetic field was set to have a uniform mass-to-flux ratio of 8 in the $x$ direction. Finally, the turbulence had been initialised with a velocity field normalised to have a Mach number of 6.7 in the central region, compatible with what is observed in Galactic star-forming clouds (Ma et al. 2022). A Kolmogorov power spectrum is then imposed to roughly reproduce the behaviour of the low-density gas filling the edges of the box, where the high temperatures make the flow sub/transonic and approximately incompressible (e.g. Porter et al. 1994). An alternative could have been to favour the low-temperature gas at the centre of the domain choosing a Burgers power spectrum, which better suits its high-Mach regime (Frisch & Bec 2000). In Sect. 4.2 we qualitatively discuss how this choice would not change the results of this work. However, we highlight that, given the idealised setup, the specific choice of the power spectrum has a marginal impact on the accuracy of the initial conditions. In effect, the mixture of high- and low-temperature gas is such that a simple power-law power spectrum cannot reproduce the turbulence spectrum characteristic of the medium<fnref n="1" />. Since the turbulent velocity field in our simulation is only an initial condition

<figure><img src="image_2.png" /><figcaption>Fig. 1. Density projections along the z-axis for two snapshots of simulation HIIR+PSJ, at 2.3 Myr (left) and 3.15 Myr (right). White dots mark the sink particles’ positions. The left panel captures the formation of the central hub, with three filaments converging into it. In the right panel we can see how the giant H II region, originating from the lower-left filament, restructured the cloud after about 0.4 Myr.</figcaption></figure>

and not continuously driven, the velocity power spectrum will naturally evolve to adopt the energy cascade that is physically consistent with the system.

To conclude, we can summarise the initial parameters in terms of the characteristic free-fall timescales of the central region, $t_{ff} \approx 1.5$ Myr<fnref n="2" />. With these initial conditions, the sound-crossing timescale of the inner region inside $r_0$ is defined by $t_{sc} = 0.15$, while the Alfvén-crossing time by $t_{ac} = 0.2 t_{ff}$. The turbulent-crossing time was set equal to the free-fall time.

## 2.3. Evolution

Given the initial conditions, which initialised the isolated cloud with a central density enhancement, the system naturally developed a hub system within a couple of megayears, with three filaments converging onto it (cf. Fig. 1, left). The four simulations evolved identically until the formation of the first star at 1.45 Myr. Afterwards, jets started to reshape the cloud and runs PSJ and HIIR+PSJ departed from the ones without jets. Lacking any feedback initially, the HIIR and NF clouds rapidly increased the total SFR – measured as the total amount of gas converted into sinks throughout two subsequent snapshots, and divided by their temporal separation (generally $\approx 0.04$ Myr). In those simulations the SFR reached its maximum at $\approx 2$ Myr (Fig. 2), while in the other two, the presence of protostellar jets lowers the SFR by a factor of 4 initially<fnref n="3" />. Between 1.8 Myr and 2.15 Myr, a few ‘low’ mass ionising stars (between 8 and $18 M_\odot$) formed in HIIR, but their associated H II regions only mildly affected the global evolution. By 2.15 Myr, while the gravity-only run settled down to a roughly constant $M_*^{tot}$, the onset of a powerful H II region, generated by a $100 M_\odot$ star, caused the SFR to drop at the same level of the other two feedback simulations. In a little more than 1 Myr the cloud was wiped out and the SFR flattened at low values. In HIIR+PSJ, the same ionising stars formed at $\approx 2.75$ Myr.

<footnote n="1">¹ It is worth noting that modern simulations can overcome this limitation by using zoom-in techniques, modelling larger galactic environments and then narrowing the focus to specific regions, using them as the initial setup for more detailed simulations (Dobbs 2015; Seifried et al. 2017). While yielding greater reliability, this approach sacrifices</footnote>

<footnote>generalisability, as the results become directly linked to the unique properties of the selected region.</footnote>

<footnote n="2">² Obtained using a mean molecular weight of 1.4, as required by the parametric cooling function.</footnote>

<footnote n="3">³ A similar impact of jets on the total SFR has also been observed in similar studies, such as Federrath (2015) – who reports a reduction of roughly a factor of 3 – Appel et al. (2023) – roughly 2 – and Murray et al. (2018) – also slightly more than 2. Although the ratio between the SFRs in NF and PSJ simulations is initially higher, agreement with these studies is recovered once both simulations reach the steady-state.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Temporal evolution of the SFR in the different simulations. The black line shows the simulation without any feedback. The run with only protostellar jets and the one with only ionising radiation are represented in blue and orange, respectively. Finally, the red line individuates the simulation with both feedback mechanisms active. The stars mark the formation of the giant H II region. The vertical dashed lines highlight the approximate time of the snapshots used in Figs. 4-5.</figcaption></figure>

Afterwards, the behaviour mirrored that of radiation-only simulation. In the right graph of Fig. 1 we can see the H II region in the process of dispersing the gas, although some triggered overdensity regions clearly arose from the shell compression. Instead, the SFR of PSJ, as that of NF, remained roughly constant up to the end of the simulation.

Therefore, it is clear that the presence of ionising stars overall reduces $\dot{M}_{\star}^{tot}$. Nevertheless, in literature, it is still greatly debated whether or not high-energy radiation can enhance the SFR and modify the star formation properties on a local scale (Pomarès et al. 2009; Elmegreen 2011; Dale et al. 2013; Roccatagliata et al. 2013; Menon et al. 2020; Wall et al. 2020). In this paper, we approached this issue from an observational perspective, investigating the underlying physical mechanisms and identifying potential biases arising from the absence of a third dimension and temporal evolution.

# 3. Methods

## 3.1. Extraction of quantities

To extract the relation between the SFR and the gas density we first projected the simulation's snapshots along the three axes at the highest resolution, giving maps of 4096² pixels. Then, we divided the column density map into a set of 500 contours, evenly log-spaced between $80 M_{\odot} pc^{-2}$ and $\Sigma_{gas} = 3 \times 10^4 M_{\odot} pc^{-2}$. As in some recent observational works (Pokhrel et al. 2021; Bieging & Kong 2022; Hu et al. 2022), the properties of a given contour level (the average density $\Sigma_{avg}$ and its area A) were obtained integrating over all the pixels above the corresponding threshold $\Sigma_{thr}$. Figure 3 gives an idea of the contours' shape, highlighting three levels at 500, 2500, 5000 $M_{\odot} pc^{-2}$ in simulation HIIR+PSJ at 2.7 Myr. We did not make distinctions between detached regions that may have formed while increasing $\Sigma_{thr}$, but we averaged all the quantities on the contour total area. Following Bieging & Kong (2022), we refer to this approach as 'cumulative' to avoid any confusion with the 'differential' one, where quantities are extracted from the region included in between two subsequent contours. We retained only contours

<figure><img src="image_4.png" /><figcaption>Fig. 3. Zoomed-in view of the central cluster of the HIIR+PSJ simulation, at 2.7 Myr. Light green, dark green, and black lines mark the contours at 500, 2500, and 5000 $M_{\odot} pc^{-2}$, respectively. The region covers an area of 1.8 pc × 1.1 pc.</figcaption></figure>

with a total surface area greater than 500 pixels, to exclude possible resolution effects on the data. As shown in Sect. 5.2 and Appendix B, this criterion was conservative in most of the cases.

Finally, we computed the SFR density $\dot{\Sigma}_{\star}$ associated with each contour as

$$ \dot{\Sigma}_{\star} = \frac{\dot{M}_{\star}}{A}, \qquad (3) $$

where $\dot{M}_{\star}$ is the total mass accreted per unit time by the sinks enclosed in the contour, calculated with the same method described in the previous section. We remark that observational studies do not have direct access to such a quantity. Indeed, the limited knowledge of YSOs' evolution leads to significant uncertainties both on their age and mass. The common approach is to assume an average mass and age for all of them, so that the SFR is proportional to the number of YSOs counted within a contour (e.g. Lombardi et al. 2014; Zari et al. 2016; Lada et al. 2017; Pokhrel et al. 2021). However, adopting this approach in our analysis would probably be affected by numerical effects. Indeed, despite the high resolution, it is not possible to carefully follow the formation of single stars and therefore to reproduce the correct initial mass function (IMF), because this would require resolution of the order of tens of AU (Lee & Hennebelle 2018). Presently, no simulation has been capable of covering such a wide range of scales while including all the physical processes treated here. We further address this issue in Sect. 5.1.

## 3.2. Free-fall timescale

The model introduced by Krumholz & McKee (2005) has the advantage of providing a physical interpretation of the law, obtaining an exponent of 1.5, similar to what is observed at galactic scales (in 2D $t_{ff} \sim \Sigma_{gas}^{-0.5}$). However, the drawback is the introduction of the new variable $t_{ff}$, particularly difficult to obtain from cloud-scale observations. The common approach is to evaluate $t_{ff}$ from the projected quantities under uniform sphere approximation (Murray 2011; Pokhrel et al. 2021; Bieging & Kong 2022). Since our purpose is to supply a possible numerical counterpart to observational studies of the KS law, we employed here the same technique. In particular, for each contour, we extracted its typical free-fall timescale estimating its volume density as

$$ \rho_{\text{avg}} = \frac{3 \sqrt{\pi}}{4} \Sigma_{\text{avg}} A^{-0.5}, \qquad (4) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the star formation–gas density relation between the different simulations, shown at four different times (*t* = 2.1, 2.7, 3.0, and 3.5 Myr), which illustrates the evolution of the cloud. In HIIR the giant H II region develops after 2.2 Myr, while for HIIR+PSJ it starts to restructure the environment at later times (> 2.75 Myr). The colour code is the same as in Fig. 2.</figcaption></figure>

and thus computed $t_{\text{ff}} = \sqrt{3\pi/32G\rho_{\text{avg}}}$, with *G* the gravitational constant.

Studies conducted on galactic surveys frequently report low values of the $\epsilon_{\text{ff}}$, ranging around 0.01 (Leroy et al. 2008, 2017; Utomo et al. 2018). However, recent works dealing with Milky Way clouds, instead, report higher values, but with a significant dispersion, going from $\epsilon_{\text{ff}} < 10^{-3}$ to $\epsilon_{\text{ff}} > 0.1$, with peaks of 0.3 (Murray 2011; Lee et al. 2016).

As mentioned in Sect. 1, more evolved models were devel-
oped in the past decade. Computing a typical timescale from
the local properties of the region, theories using the multi-
free-fall time proved to provide valid results on galactic and
extragalactic scales. However, we highlight that our intent is to
compare the star formation properties of the inner structures of
molecular clouds as recovered from numerical simulations and
observations. In this context, we are not aiming to utilise the
most advanced model to calculate $\epsilon_{\text{ff}}$, but rather follow the very
same approach used in the works based on observation we are
comparing to.

4. Resulting star formation relations

4.1. Kennicutt–Schmidt relation

We show the resulting KS relation obtained for our clouds
in Fig. 4. To facilitate the comparison between simulations,
we scaled the y-axis by $\Sigma_{\text{avg}}^2$, which represents the character-
istic slope observed in the radiation simulations. This is also
consistent with typical values reported in observational stud-
ies conducted on cloud scales (see e.g. Gutermuth et al. 2011;
Lada et al. 2013; Pokhrel et al. 2021). The plots show the projec-
tion along the z-axis, but comparisons with other lines of sight
yielded similar conclusions (see our Appendix A and Khullar

et al. 2019). We selected four specific snapshots to represent key
stages in the simulations’ evolution:

- At 2.1 Myr: by this time the SFR of run NF has almost reached the steady state (Fig. 2), while cloud HIIR has not been affected by the eruption of the large H II region.

- At 2.7 Myr: here, due to the presence of the 100 $M_{\odot}$ star, the SFR of the radiation-only simulation decreased and reached the level of the runs with jets. The same star has yet to be formed in HIIR+PSJ.

- At 3.0 Myr: at this time, the four clouds all present signif-
icant differences between each other, and the main H II region
has started to restructure the full-feedback run (HIIR+PSJ).

- Finally, we plot the situation at the very last moments of
the radiation simulations, at 3.5 Myr.

The plots reveal that the different feedback prescriptions lead to significant differences in the clouds’ evolution. In the early stages, when only a few small H II regions developed, jets dominate as the primary feedback mechanism. Therefore, at 2.1 Myr, we observe strong similarities between the runs NF and HIIR (without jets) and PSJ and HIIR+PSJ (with jets). Similarly to what is displayed in Fig. 2, we see that stellar ejecta reduce the global SFR by a factor of $4^{4,5}$, without sensibly modifying the shape of the relation.

By 2.7 Myr the situation has changed substantially. Overall,
the KS relation has become an almost pure power law in all
simulations. Towards the low-density edge, a clear change of

<footnote n="4">⁴ The leftmost side corresponds to the widest contours, which include all the others. From those we can recover the average behaviour of the cloud.</footnote>

<footnote n="5">⁵ Notice that this factor is greater than the loss of accretion due to the expulsion of material, which alone contributes to decreasing the total SFR by a factor of 1.5.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the modified KS relation. The colour code is the same as in Fig. 2. In the graph at 3 Myr, we report the values of $\epsilon_{ff}$ of the 12 Galactic clouds studied in Pokhrel et al. (2021), excluding contours with $\Sigma_{avg}$ above their reliability criterion (grey lines). The grey shaded areas represent the uncertainties they gave for the $\epsilon_{ff}$ measures. The axis ratios are kept the same as in Fig. 4 to facilitate a visual comparison of slopes with the respective counterparts in the KS relation.</figcaption></figure>

slope is visible in runs with feedback, while it is almost invisible in run NF. Probably, this broken power law identifies the point at which stellar feedback counteracts the global collapse, and the cloud's KS relation departs from that typical of a free-fall motion. Moreover, we can note HIIR+PSJ showing a slightly steeper relation compared to the jets-only simulation. Interpreting the graph, this means that the relative fraction of gas mass that is converted into stars per unit time increases at high densities when ionising radiation is included, even though at this time the total SFR is similar (Fig. 2). This behaviour is even more accentuated comparing runs HIIR and NF, where the first achieves higher values of $\Sigma_{\star}$ at high densities, with an exponent of its KS relation of ≈2.

At 3.0 Myr, the different exponent of the HIIR+PSJ and PSJ KS relation became even more evident. Another visible feature is the departure from the power-law behaviour of run HIIR at $\Sigma_{avg} \approx 10^3 M_{\odot} pc^{-2}$. Given the sudden decline, we can deduce that the main H II region, dispersing the gas, sensibly lowered the column density around a sink formation site. Then, when $\Sigma_{thr}$ exceeds the maximum column density of that region, the observed SFR suddenly drops. Run NF, instead, exhibited an increase in $\Sigma_{\star}$ towards higher column densities. In the simulation, this increase coincides with the moment in which a clump, formed in one of the filaments, converges into the central hub.

Finally, at 3.5 Myr, the large H II region has fully developed in the simulation HIIR+PSJ and its $\Sigma_{\star}$ versus $\Sigma_{avg}$ relation mirrors that of the radiation-only simulation. Its slope, as well as the one from HIIR, closely approaches $\Sigma_{avg}^2$. With the system reaching a steady state, the NF run does not display

significant differences compared to the previous graph, while PSJ experiences a further steepening of the relation.

## 4.2. Modified Kennicutt–Schmidt relation

The evolution of $\epsilon_{ff}$ across time (Fig. 5), obtained inverting Eq. (1), closely reproduces that seen in Fig. 4. However, the differences in both slope and scatter between the PSJ and radiative simulations are less pronounced. In effect, one of the major strengths of the modified relation is its reduced scatter among different objects, compared to the traditional KS law (Hu et al. 2022). Observational datasets typically consist of multiple clouds at various evolutionary stages, characterised by different strengths of feedback mechanisms within the complex. Given that our simulated clouds feature different feedback prescriptions, the smaller dispersion observed in the four $\epsilon_{ff}$ values aligns well with expectations based on observational results.

As in Fig. 4, the curves from radiative simulations tend to achieve higher values towards the high-density end of the graphs. Additionally, the $\epsilon_{ff}$ of these simulations appears almost flat, while in runs NF and PSJ it shows a steady decline of $\epsilon_{ff}$ as the density increases. Indeed, this qualitative result is confirmed when fitting a power law such as

$$ \epsilon_{ff} \sim \Sigma_{avg}^{\gamma} \qquad (5) $$

At each snapshot, we fitted<fnref n="6" /> the relation using only contours that satisfy the 500 pixels contours and have $\Sigma_{avg} > 300 M_{\odot} pc^{-2}$.

<footnote n="6">⁶ We employed the Python built-in function `curve_fit`.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the slope of $\epsilon_{\text{ff}}$ across the simulations. The shaded areas outline the $1\sigma$ standard deviation of the fitted exponent as obtained from the three different axes of projection. The stars mark the formation of the main H II region, as in Fig. 2</figcaption></figure>

This allowed us to exclude the different relation visible at lower densities (i.e. larger scales), which might be affected by the initial condition of the cloud. Indeed, the evolutionary timescales of these regions are much longer, and the system employs more time to relax from the initial conditions here. We repeated this procedure for all three projection axes and assigned to the average $\gamma$ an uncertainty given by the standard deviation from the three lines of sight (the errors from the fit procedure were completely negligible).

Figure 6 shows the evolution of the exponent $\gamma$ across the simulation. Again, we see that only simulations with radiation ultimately set on $\gamma$ values close to zero. A $\epsilon_{\text{ff}}$ independent of $\Sigma_{\text{avg}}$ has also been reported in nearby Galactic cloud by Pokhrel et al. (2021), although in their dataset the value of $\gamma$ scatters between -0.33 (AFGL 490) and 0.06 (Cygnus-X)<fnref n="7" />.

Nevertheless, compared to their work, we found sensibly higher values of $\epsilon_{\text{ff}}$. For comparison, we overplotted the curves they obtained for the 12 Galactic clouds in the graph at 3 Myr in Fig. 5, and we can see how the two measures differ by a factor of ≈ 5. While the influence of the SFR prescription cannot be ruled out entirely, it is improbable that such a difference is attributable entirely to the code setup<fnref n="8" />. Indeed, we expect the total simulated SFR to be reasonably accurate in runs with jets<fnref n="9" />, and the sink spatial distribution reliable at least up to $\Sigma_{\text{thr}} > 3000 M_{\odot} \text{ pc}^{-2}$ (roughly 3.5 in the graph, cf. with Fig. 8). Moreover, even though our simulations could not accurately follow the dynamical evolution of sinks due to the softening of gravity, Lombardi et al. (2013) and Lada et al. (2013) show that in real clouds the diffusion process due to N-body interactions does not contribute significantly to the observed $\Sigma_{\star}$ at low densities. Finally, we outline that adopting a Burgers power spectrum as initial condition for turbulence would probably not solve the tension observed with the $\epsilon_{\text{ff}}$ reported in Pokhrel et al. (2021). Indeed, a steeper power spectrum means that much less kinetic energy is stored

in small scales, leading to more coherent large-scale flows. Therefore, imposing a Burgers spectrum in the velocity field would have enhanced the formation of bound structures (Clark & Bonnell 2005), likely increasing the cloud's SFR and consequently the discrepancy with the SFRs obtained from the observational data.

Therefore, we looked at three potential biases that may influence observational measurements. We emphasise right away that our intent is to understand how different biases could interact with the feedback mechanisms to affect the results. Although full agreement with Pokhrel et al. (2021) $\epsilon_{\text{ff}}$ measures is recovered, giving a precise quantitative prediction of the impact of these biases is beyond the purpose of this paper.

# 5. Impact of biases on the observed relation

We examined, in our simulations, the effect of three biases that affect observations: the assumption that the SFR is proportional to the number of Class 0/I YSOs detected, the impact of resolution on column density maps, and the presence of embedded stars not visible to the telescope due to extinction. In this section we aim to understand their impact on the observed $\epsilon_{\text{ff}}-\Sigma_{\text{avg}}$ relation. In Fig. 7, we plotted the effects of these biases separately, using the snapshot at 3 Myr as a reference (blended curves in the background). At that time, the contribution of ionising radiation in HIIR+PSJ was not as dominant as in HIIR yet. This allowed us to study how these biases intertwine with the different levels of feedback.

## 5.1. YSO counts

The first bias we considered is caused by the uncertainty on detected YSOs' age and mass. Since observations cannot access the temporal evolution of the star-forming clouds, observers rely on theoretical models to constrain these properties in young stars. As mentioned in Sect. 3, the usual procedure is to assign a mean mass and age to each of them. The mass corresponds to the average mass of the IMF, which is typically around $0.5 M_{\odot}$. For Class 0/I objects, the age is commonly assumed to be 0.5 Myr (Dunham et al. 2015; Megeath et al. 2022). With these values, the SFR of each contour is $\dot{M}_{\star} = 1 M_{\odot} \text{Myr}^{-1} \times N_{\text{YSO}}$, where $N_{\text{YSO}}$ is the number of young stars it contains. Imitating this technique, we counted only the sinks younger than 0.5 Myr. However, we explained in Sect. 3 that the IMF in our cloud is not well sampled due to insufficient resolution. Indeed, the average stellar mass resulted to be higher than the expected $0.5 M_{\odot}$, with $\langle m_{\text{sink}} \rangle \approx 0.75 M_{\odot}$ in the feedback runs and $2 M_{\odot}$ in NF at the end of the simulation. Using a coefficient of $1 M_{\odot} \text{Myr}^{-1}$ in the above equation would underestimate the total SFR by a factor of $\langle m_{\text{sink}} \rangle/0.5 M_{\odot}$. Therefore, we corrected this by assigning to each contour an internal SFR given by $\dot{M}_{\star} = \langle m_{\text{sink}} \rangle/(0.5 \text{ Myr}) \times N_{\text{YSO}}$.

The first plot in Fig. 7 (upper left) demonstrates that this particular approximation well reproduces the true SFR of the cloud. Except for NF, the detected $\epsilon_{\text{ff}}(\Sigma_{\text{avg}})$ preserves good agreement at all contours. Minor deviations are observed at high-density contours, probably due to random fluctuations in the instantaneous SFR. In these regions, $t_{\text{ff}}$ becomes much smaller than the assumed build-up time of 0.5 Myr, so that the integration time becomes much longer than the dynamical timescale (Pokhrel et al. 2021).

<footnote n="7">⁷ This can be recovered from the second model fitted in Table 1 (sixth column), since $\gamma = a - 1$.</footnote>

<footnote n="8">⁸ Our findings are not the first to reveal a mismatch between the observed and simulated SFRs and $\epsilon_{\text{ff}}$. Specifically, Federrath (2015) has also highlighted that numerical simulations often result in higher SFRs than those reported in observational studies.</footnote>

<footnote n="9">⁹ The lack of gas expulsion from sinks in HIIR and NF causes an overestimate of the SFR by a factor of 1.5.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Effect of four different observational biases on the $\epsilon_{ff} - \Sigma_{avg}$ relation of the clouds at 3 Myr. The blended curves are the same as in the 3 Myr panel of Fig. 5, and the grey lines are the 12 Galactic clouds studied in Pokhrel et al. (2021).</figcaption></figure>

## 5.2. Spatial resolution

To reproduce the second bias, we smoothed the column density maps, lowering the resolution of our ‘observation’. This bias cannot impact the measured SFR, but rather affects the observed cloud structure<fnref n="10" />. Placing the cloud at a distance<fnref n="11" /> of 650 pc, we reproduced the physical resolution of the Herschel satellite in the 500 µm band (36 arcsec), convolving the column density maps with a Gaussian beam with a standard deviation of 0.05 pc (equivalent to a full width half maximum of 0.11 pc). We did not modify the criterion given in Sect. 3, because now the threshold does not have to define a limit for the physical reliability of the data. Moreover, in this way, we could analyse how a low spatial resolution can affect the measured $\epsilon_{ff}$ if this limitation is not carefully considered.

Smoothing the maps produces what is visible in the upper right graph of Fig. 7. Except HIIR, the curves of the simulations almost overlap with the original ones below $1000 M_{\odot} pc^{-2}$ (presented in Fig. 5). Then, an increase is visible at the right end in HIIR+PSJ and PSJ, while this is evident at low $\Sigma_{avg}$ in the radiation-only run. This behaviour can be easily explained by examining how the function $A(\Sigma_{avg})$ enters the relation.

Combining Eqs. (3), (1), and (4), we can express $\epsilon_{ff}$ as<fnref n="12" />

$$ \epsilon_{ff} = \frac{\dot{\Sigma}_{\star} t_{ff}}{\Sigma_{\text{avg}}} \sim \frac{\dot{M}_{\star}}{A} \cdot \Sigma_{\text{avg}}^{-3/2} A^{1/4} \sim \dot{M}_{\star} \Sigma_{\text{avg}}^{-3/2} A^{-3/4}. \quad (6) $$

It follows that the cloud structure – $A(\Sigma_{avg})$ – plays a central role in defining the behaviour of $\epsilon_{ff}$ across the column density contours. As we approach the resolution limit, the integrated probability distribution function of $\Sigma_{avg}$ (usually referred to as integrated N-PDF) drops sharply, leading to an increase in the efficiency per free-fall time thanks to the negative exponent of A. This causes the knee visible at high densities.

In effect, by plotting the contour area as a function of their $\Sigma_{avg}$ at high and low resolution in Fig. 8, we can identify the same feature. Moreover, we notice that the adopted threshold of 500 pixels (dashed line) effectively safeguards the data from resolution effects, as in the first graph the exponential decrease occurs beyond this point. Degrading the resolution moved this feature to lower $\Sigma_{avg}$, which caused the curves to bend towards the high-density end of Fig. 7.

However, this alone cannot explain the full picture. Indeed, we see a steepening of the power-law tail of integrated N-PDF even at low densities. Technically, such variation should not be expected, as the smoothing should not significantly affect

<footnote n="10">¹⁰ Resolution applied to YSOs detection constitutes instead a different issue. In high $\dot{\Sigma}_{\star}$ zones, the crowding effect could lead to consistent decreases in completeness. Megeath et al. (2016, 2022) showed how accounting for this number of YSOs could increase by 30% in the Orion complex.</footnote>

<footnote n="11">¹¹ That is the average distance of the clouds considered in Pokhrel et al. (2021) study.</footnote>

<footnote n="12">¹² Notice that the usual $t_{ff} \sim \Sigma_{avg}^{1/2}$ imply the presence of a constant scale width much smaller than the surface. While true when looking at external galaxies, this is no longer valid for clouds, where the third dimension is ‘comparable’ to the others and its dependence must be included.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Cumulative area of the contours as a function of their average column density at 3 Myr. The two images are obtained from the original column density map (top) and the one smoothed with a Gaussian beam with 0.05 pc standard deviation (bottom). As in Fig. 7, the blended lines in the bottom panel are the ones plotted at 0.0073 pc resolution. The dashed line marks the area corresponding to the 500 pixel criterion.</figcaption></figure>

areas much larger than the beam size<fnref n="13" /> (see e.g. Tassis et al. 2010; Schneider et al. 2015). Most probably, the pitfall lies in the integration over the whole contour. As explained in Sect. 3, the method we applied does not take into account separated structures within the same column density level. Consequently, resolution can have a much stronger impact on the observed N-PDF of a highly fragmented structure. At 3 Myr, the HIIR cloud has been significantly restructured by the H II region, resulting in the destruction of the central hub and the formation of several disconnected high-density clumps. The smoothing gives rise to a completely different γ evolution.

### 5.3. Instrument sensitivity

Finally, we accounted for the possibility that some YSOs may be too heavily extinct to be detected in the inner regions of the cloud. We considered the observations to be taken with the Spitzer Infrared Array Camera (IRAC) 4.5 µm filter, which is extensively adopted in SFR measurements (see e.g. Evans et al. 2009; Gutermuth et al. 2011; Lombardi et al. 2013). For each sink, we evaluated its apparent magnitude as

$$m_{4.5} = M_{\text{bol}} - BC_{4.5} + 5 \log \left( \frac{d}{10 \text{ pc}} \right) + A_{4.5}, \quad (7)$$

where $M_{\text{bol}}$ is the bolometric luminosity of the sink, $BC_{4.5}$ its bolometric correction for IRAC 4.5 µm filter, $d$ the distance of the cloud, which we set to 650 pc, and $A_{4.5}$ the extinction at 4.5 µm.

To recover the needed stellar properties such as the luminosity ($L$), radius ($R$), and bolometric correction, we compared the age and mass of each sink with the evolutionary tracks obtained from the MIST web tool<fnref n="14" />. However, considering that in Class 0/I objects the luminosity due to gas accretion, $L_{\text{acc}}$, is often a consistent fraction of their total luminosity $L_{\text{tot}}$ (Antoniucci et al. 2008; Fiorellino et al. 2021), we corrected the sink’s bolometric luminosity as $L_{\text{tot}} = L + L_{\text{acc}}$ (Fiorellino et al. 2023). Following Hartmann et al. (1998), we define $L_{\text{acc}}$ as

$$L_{\text{acc}} = 0.8 \frac{GM}{R} \dot{M}, \quad (8)$$

where $M$ and $\dot{M}$ are the sink mass and its accretion rate averaged over ≈ 0.05 Myr. Finally, we obtained the last ingredient $A_{4.5}$ evaluating the amount of gaseous material in front of the sinks, with the conversion $A_{4.5} = 0.572 \cdot A_{\text{K}}$ provided in Ascenso et al. (2013), and adopting $\Sigma_{\text{gas}} = 197 M_{\odot} \text{pc}^{-2} \cdot A_{\text{K}}$ (hence, $\Sigma_{\text{gas}} = 384 M_{\odot} \text{pc}^{-2} \cdot A_{4.5}$; Lada et al. 2013).

To be detected, we required sinks to have $m_{4.5} < 14.5$ mag, which is the sensitivity limit at 4.5 µm reported in the Spitzer/IRAC Candidate YSO (SPICY) catalogue (Kuhn et al. 2021). However, we stress once again that, with these simplified prescriptions, our goal is to understand how those biases could affect the observed $\epsilon_{\text{ff}}-\Sigma_{\text{avg}}$ relation. Indeed, observers employ different instruments at different wavelengths to limit the impact of these biases (Harvey et al. 2007; Megeath et al. 2016). These corrections are particularly demanding to reproduce in simulations and require specific treatments that are beyond the scope of our study.

The effect of the instrument sensitivity is shown in the third graph of Fig. 7 (lower left). A notable feature is the reduction of $\epsilon_{\text{ff}}$ average value. Even if caution has to be used in the interpretation, the resolution of the tension with Pokhrel et al. (2021) data could imply that the physics implemented in PSJ and HIIR+PSJ is sufficient to describe the mechanisms involved in setting the $\epsilon_{\text{ff}}$ in molecular clouds. In our simulations, adding a detection threshold has a strong impact on the observed SFR, leading to an underestimation of a factor of 3 of $\epsilon_{\text{ff}}$ in the full-feedback run, and ≈ 5 in the jet-only run. This difference is solely caused by the different distribution of YSOs, as supported by the fact that this bias has a softer impact on the HIIR simulation. Indeed, by 3 Myr the main H II had already significantly restructured the cloud, dispersing a consistent part of the gas, reducing the extinction experienced by stars.

Another consequence of this bias is the rapid fall of the relation at high $\Sigma_{\text{avg}}$. Being a log-log graph, the slope of $\epsilon_{\text{ff}}$ is directly affected by the fractional variation of $\dot{M}_{\star}$ at each contour. Most of the stars – and so the majority of the SFR – reside at high column densities. Therefore, all low-$\Sigma$ contours have comparable $\dot{M}_{\star}$. Similar drops are also observed in Pokhrel et al. (2021) clouds when $\Sigma_{\text{thr}}$ exceeds their reliability criterion (defined by the column density at which $t_{\text{ff}} > 0.5$ Myr). It is important to highlight that this reduction is specific to the cumulative method. Since $\Sigma_{\star}$ is obtained through the integration over the entire contour, missing YSO in the central regions lowers the SFR of every contour, even those with low $\Sigma_{\text{avg}}$. In contrast, a differential approach like that used in Evans et al. (2014) would not be

<footnote n="13">¹³ With little algebra, it is possible to demonstrate that the exponent of the N-PDF computed in logscale is the same of the one obtained from the integrated N-PDF (Lada et al. 2017).</footnote>

<footnote n="14">¹⁴ https://waps.cfa.harvard.edu/MIST/interp tracks.html</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 9. SFR of density contours as a function of their threshold column density at 3.0 Myr. The data are normalised to the total SFR at that time. The dotted part of the lines identifies isosurfaces with an area smaller than 500 pixels. The right panel shows the same distribution, but as it appears after the application of all three biases explained in the text – YSO counting for the SFR, resolution lowered to 0.05 pc, and the masking of extinct stars. The dashed green lines reproduce the model from Lada et al. (2017), with a threshold of 2000 M$_{\odot}$ pc$^{-2}$ and a stellar probability distribution function of PDF$_{*}$ ∝ Σ$^{0.7}$thr.</figcaption></figure>

affected, but it would require rougher assumptions for comput-
ing t$_{ff}$ and would suffer from the lower statistics of YSOs found
at low densities.

We can show this phenomenon more accurately if we extract
the relation between Ṁ<sub>⋆</sub> and the contour column density. In
Fig. 9a, we plotted the evolution of Ṁ<sub>⋆</sub> for different contour
levels at the maximum resolution, as a function of the column
density threshold of that contour Σ<sub>thr</sub>. With this change of
x-axis, we can provide a comparison with the model by Lada
et al. (2013), shown in the graph as a dashed green line.
According to this work, the SFR included in each contour can
be described as

$$
\dot{M}_{\star} = \alpha_{\mathrm{N}} \int_{\Sigma_{\mathrm{thr}}}^{\Sigma_{\mathrm{max}}} \mathrm{PDF}_{*} \, \mathrm{d}\Sigma_{\mathrm{gas}}, \qquad (9)
$$

where αN is the normalisation constant and PDF* is the proba-
bility of finding a star at a given gas column density. As derived
from the California Cloud in their work, we used

$$
\mathrm{PDF}_{*} \propto \begin{cases} \Sigma_{\mathrm{gas}}^{\phi} & \text{if } \Sigma_{\mathrm{gas}} > \Sigma_{\mathrm{min}} \\ 0 & \text{otherwise} \end{cases}, \qquad (10)
$$

with $\phi = -0.7$ and $\Sigma_{\min} = 100 M_{\odot} \text{pc}^{-2}$, and integrated up to a column density $\Sigma_{\max} = 2000 M_{\odot} \text{pc}^{-2}$ ($A_{\mathrm{K}} \approx 10 \text{ mag}$).

The net differences between our simulations and Eq. (10) are
clear and quantitative. In fact, attempting a fit leaving Σ<sub>max</sub> and
φ as free parameters and Σ<sub>min</sub> equal to the threshold at which
the first accreting star is removed, leads to φ ≈ 1 in all the sim-
ulations, while for the four clouds studied Lada et al. (2017)
this varies between -1 and -0.3. This picture changed after
the application of all the observational biases described above
(Fig. 9b). Qualitatively, the curves now resemble the model
much more closely. Repeating the fit of Eq. (10), we found that
only HIIR+PSJ and PSJ converged to reasonable values of φ
(-0.32 ± 0.02 and -0.42 ± 0.02, respectively)<sup>15</sup>.

We can evaluate the derivative of Eq. (9) to understand how
the sensitivity changes the contribution of Ṁ*(Σ<sub>thr</sub>) ~ Ṁ*(Σ<sub>avg</sub>)

to the observed KS relation. Assuming that Eq. (10) holds at
every density contour, with little algebra we arrive at

$$
\frac{d \ln \dot{M}_{\star}}{d \ln \Sigma_{\text{thr}}} (\Sigma_{\text{thr}}) = -(\phi + 1) \frac{\chi^{\phi+1}}{1 - \chi^{\phi+1}}, \qquad (11)
$$

with χ = Σthr/Σmax. At lower densities, two ingredients con-
tribute to the initial slope of Ṁ⋆. The first is exponent φ, which
Lombardi et al. (2013) show to be the sum between the exponent
−q of the A(Σavg) function and the exponent of the KS relation
as obtained from the differential approach. The steeper this last
relation, the more stars will be present at higher densities and
hence the shallower the logarithmic derivative will be. The sec-
ond is χ, which is the combination of the scale at which Eq. (10)
ceases to be valid, and the sensitivity limit of the instrument.

5.4. Biases combined

To conclude, in the last panel of Fig. 7, we plotted ϵ<sub>ff</sub> as seen
after the application of all these biases on our column density
maps. First, we masked sinks too extinct to be detected. Then,
we smoothed the column density maps. Finally, we repeated the
same analysis described in Sect. 3, using the new YSO count
prescription to evaluate ȳ<sub>⋆</sub>.

The resulting relations are similar to what we would obtain
from the simple sum of all the effects, meaning that any coupling
between the biases has little effect on the resulting relation. How-
ever, a notable exception arises when looking at the curve of NF.
In this panel, ϵ<sub>ff</sub> is much higher than the one observed in the
‘Sensitivity’ graph. This cannot be due to the smoothing of the
column density map, as we said that this does not affect the total
SFR detected. Therefore, the YSO counting must increase the
observed SFR, although the first panel clearly shows that this
decreases the average ϵ<sub>ff</sub> when taken alone.

The answer resides in the same argument we used in
Sect. 5.1. Counting YSOs accurately reproduces the SFR only
when the correct average mass is used. However, if the mass dis-
tribution varies as we approach the densest regions of the cloud,
then a new bias is introduced. In effect, as we increase Σ<sub>thr</sub> we
observe a consistent increase in the average stellar mass of the
contours. In NF, the average sink mass in the hub increases sig-
nificantly, and thus ⟨m<sub>⋆</sub>⟩ overestimates the mean mass of YSOs

<footnote n="15"><sup>15</sup> Run NF, instead, returns a value of φ = 0.68 ± 0.03, while in run HIIR the fit does not converge.</footnote>
found at low densities while underestimating the ones found at the highest ones. Hence, if we remove the stars in the centre of the cloud, this approximation tends to overestimate the SFR observed.

We remark that in our cloud, this phenomenon is probably attributable to numerical effects (see Sect. 3), since in all the simulations containing feedback $\langle m_{\star} \rangle$ remains fairly constant even at high $\Sigma_{thr}$. Nevertheless, this offers insights into the interpretation of observations as well. Whether the mass function varies as a function of the local density is a subject still under debate (Raboud & Mermilliod 1998; Stolte et al. 2002; Bonnell 2008). Moreover, in young massive clusters, the mass segregation process takes place in very short timescales, possibly leading to an effect of this kind (Portegies Zwart et al. 2010; Allison et al. 2010; Suin et al. 2022).

# 6. The shape of the relation

## 6.1. Cloud structure

In the previous section, we have found clues that the function $A(\Sigma_{avg})$ plays a central role in shaping the relation. Figure 8 demonstrated that the effects produced by low resolution on the observed cloud structure directly translate into similar features in the $\epsilon_{ff}$ versus $\Sigma_{avg}$ graph. Simultaneously, Fig. 9 shows that the more we can access the SFR in the highest densities of the cloud, the more $\dot{M}_{\star}$ appears flat at lower densities. This is because most of the SFR happens in the densest regions. Consequently, the shape of the relation at low-mid densities is mostly governed by the slope of the $A(\Sigma_{avg})$, while at high densities it intertwines with the loss of stars between contours. Which of the two dominates in this regime may vary from cloud to cloud (e.g. in the last panel of Fig. 7, compare HIIR with the other runs).

In Sect. 5.2, we find that the presence of ionising feedback does play a role in restructuring the cloud, changing the observed area-to-column-density relation. Even at the maximum resolution, in Fig. 8, some differences are evident between the various setups. In the radiation runs, the compression caused by the expansion of the H II region fronts is not sufficient to overcome the gas dispersion, and comprehensively the total amount of dense gas decreases, leading to a steeper relation. To measure this effect quantitatively, we fitted a power-law relation

$$A \sim \Sigma_{\text{avg}}^{-q} \qquad (12)$$

at each snapshot. We performed the fit with the Python built-in function `curve_fit`, using only the contours having $\Sigma_{avg} > 400 M_{\odot} pc^{-2}$ (log $\Sigma_{avg} \approx 2.6$) and an area greater than 500 pixels. The first constraint allowed us to exclude from the fit the largest scales of the simulations, which are dominated by the global collapse and thus display an exponent $q$ similar to that of NF. In Fig. 10, we plotted the fit results, averaged over the three axes of projection. The uncertainties associated with the fit are negligible, so the shaded area represents the dispersion between the three lines of sight. Without any support from protostellar jets, the runs NF and HIIR collapsed faster and achieved lower values of $q$ (more dense gas at high $\Sigma_{avg}$). The onset of large H II regions froze the exponent (at $\approx 2.3$ Myr in HIIR and $\approx 2.8$ Myr in HIIR+PSJ), which settled to a value of roughly 2 (>1.85 in HIIR+PSJ, $\approx 2.05$ in HIIR).

Conversely, in NF the exponent continued to decrease down to roughly 1.6, close to what is reported and expected in gravoturbulent simulations without feedback (Burkhart et al. 2015, 2017; Chen et al. 2018). Notably, the run with only jets achieves a shallower relation than the run with only gravity, reaching values as

<figure><img src="image_11.png" /><figcaption>Fig. 10. Temporal evolution of the fitted exponent $q$, averaged over the three axes, in the four simulations. The colour code is the same as in Fig. 2. The shaded areas represent the standard deviation of the sample due to the differences between the lines of sight. The uncertainties on the single fits are negligible. As in Fig. 2, the stars highlight the moment of the eruption of the main H II region.</figcaption></figure>

low as 1.4. We notice that the values of $q$ between NF and PSJ cross almost exactly at 3 Myr, so that the steeper relation is not visible in Fig. 8. However, looking at Fig. 2, we see that this excess of dense gas is not rewarded by an increase in the total SFR. This is probably because compressed material is ejected at high relative velocities relative to the surrounding gas flow.

We can compare these findings with what is usually observed in real molecular clouds<fnref n="16" /> (cf. Lombardi et al. 2015; Spilker et al. 2021). Our simulated clouds attained slopes similar to the observed ones, although it appears that values of $q$ near or above 2 are more frequent in our Galaxy. Moreover, the fact that only clouds with radiation attained a similar exponent may represent a central result of this work. Indeed, Abreu-Vicente et al. (2015) conducted an extensive study on the N-PDF distribution in 195 Galactic molecular clouds, concluding that clouds with H II regions naturally develop a power-law tail with an index -2. The authors suggest this slope to be caused by the formation of an isothermal-sphere-like structure. Indeed, based on the work from Kritsuk et al. (2011), they show that for a spherical density distribution of the form $\rho_{gas} \sim r^n$, with $r$ the distance from the centre, a relation between $n$ and $q$ can be established, such that $n = 1 + 2/q$. Setting $q = 2$, we recovered the isothermal sphere distribution. Looking at Fig. 8, the presence of a mechanism that pushes $q$ towards this value is also suggested by the smaller dispersion between the line of sights that characterises the radiative runs (although a larger sample of simulation would be needed to confirm this result statistically). This could potentially mean that radiation feedback plays a fundamental role in setting the density structure of clouds. This could provide valuable insights to distinguish between different classes of cloud.

Interestingly, the same study from Abreu-Vicente et al. (2015) detects that on average the star-forming clouds observed without H II regions display slopes steeper than the one reported here. This seems to contrast with our results. Nonetheless, such a steep slope is visible at the early stages of the simulations (left-hand side of Fig. 8). Thus, a possible explanation is that the clouds that have such a high $q$ are simply younger than the others. Possibly, those clouds will proceed to evolve a shallower slope,

<footnote n="16">¹⁶ Notice that Spilker et al. (2021) fits the exponents $\alpha$ of the N-PDF in logspace. Hence, their exponent $\alpha$ is identical to our $q$ (footnote 13).</footnote>
but eventually form a massive star that interrupts the process as happened in HIIR+PSJ.

## 6.2. Impact of feedback

Figure 10 gives us the necessary instruments to interpret the behaviour of the observed KS relations in Figs. 4–5, since we found in the previous section that $\dot{M}_{\star}$ has very little impact at low column densities. In the two figures, the radiative simulations show a slope nicely falling onto a $\dot{\Sigma}_{\star} \sim \Sigma_{\text{avg}}^2$ in the first and $\epsilon_{\text{ff}} \sim \text{const}$ in the second. This behaviour can be perfectly explained by the reported $q \approx 2$. Concerning the KS relation, combining Eq. (3) with the $A(\Sigma_{\text{avg}})$ in Eq. (12), it is clear that a nearly constant $\dot{M}_{\star}$ leads to $\dot{\Sigma}_{\star} \sim \Sigma_{\text{avg}}^q$. Doing the same for Eq. (6), it is easy to find that under the same $\dot{M}_{\star} \sim \text{const}$ approximation, a value of $q \approx 2$ leads to a constant $\epsilon_{\text{ff}}$. Looking closely we can find that Figs. 6–10 closely resemble each other. At the same time, the weaker dependence on $A$ of Eq. (6) allows the relation to be more resilient to the effect of feedback, explaining the smaller deviations between the simulations compared to the differences found in the KS relations (see Sect. 6.1).

This also suits the shallow behaviour of HIIR and HIIR+PSJ in the 3D analysis (Fig. B.1). With little algebra, it is possible to show that if the same assumption of spherical symmetry is valid, then $V \sim \rho_{\text{gas}}^{-3/q}$. If $q \approx 2$, the influence of the volume in the tri-dimensional version of Eq. (6) is expected to be negligible, since

$$ \epsilon_{\text{ff, 3D}} = \frac{\dot{M}_{\star}}{V\rho_{\text{gas}}} t_{\text{ff}} \sim \dot{M}_{\star} \rho^{\frac{3q}{2+q}-\frac{3}{2}}, \quad (13) $$

with the exponent going to zero for $q=2$.

Summing up all these findings, we conclude that when ionising radiation is included in the simulation, the cloud systematically develops a KS relation more weighted towards higher column densities, as well as a higher $\epsilon_{\text{ff}}$ and $\epsilon_{\text{ff, 3D}}$ – although with a reduced impact. Thus, although radiation negatively impacts the SFR on the largest scale and reduces the amount of dense gas available, the remaining dense gas manifests a higher $\epsilon_{\text{ff}}$.

Unfortunately, the interpretation of the observed trend is not unique. One possible explanation is that the expansion of ionised gas promotes star formation at high densities in an active way. In this scenario, the radiation front would create favourable conditions for the development of new star formation sites, compressing the surrounding gas until it becomes gravitationally unstable (collect and collapse mechanism; Elmegreen & Lada 1977), creating converging shocks at clump surfaces capable of enhancing significantly their core density (radiation driven implosion; Sandford et al. 1982), or acting as a compressive driver for local turbulence, (as suggested by Menon et al. 2020, 2021). Alternatively, it could be suggested that massive stars wipe out only the gas that is not sufficiently gravitationally bound, leaving only the dense gas that would have formed stars regardless. This again would lead to deriving a higher star formation efficiency in the region, since the local SFR would roughly remain the same, but the observed amount of dense gas would be lower. Which of the two processes prevails is unclear and subject to intense studies in the field. Either way, the outcome of this interplay produces observable features in our cloud, which could aid in interpreting observational data. For instance, in a sample of real clouds, the differences in their evolutionary stage or variations in the feedback mechanisms' strengths could contribute to the observed scatter in the derived KS laws.

Indeed, Figs. 4–5 show that the relations undergo significant variations over a timescale of a few million years, consistently

<figure><img src="image_12.png" /><figcaption>Fig. 11. Time average of $\log \epsilon_{\text{ff}}$ at 250, 500, 1000, 2500, and $5000 M_{\odot} \text{pc}^{-2}$ in the various simulations. The data have been slightly displaced to make the comparison clearer. The different markers highlight different axes of projections, and the error bars show the dispersion of $\epsilon_{\text{ff}}$ measured at a given column density due to its evolution across the simulation. The effects from the line of sight and stellar feedback are not enough to overcome this dispersion, such that in a set of clouds at different evolutionary stages they will not be visible if the dependence from $\Sigma_{\text{avg}}$ is not taken into account.</figcaption></figure>

with what is reported in Grudić et al. (2018). As time progresses, the four extracted relations intersect at different column densities. This has implications for observational studies that extract single measurements of $\epsilon_{\text{ff}}$ or $\dot{\Sigma}_{\star}$ from each cloud. These analyses often adopt a common minimum $\Sigma_{\text{thr}}$ to define the extent of the clouds. Our results emphasise that achieving similar values of $\epsilon_{\text{ff}}$ or $\dot{\Sigma}_{\star}$ at a specific column density contour does not necessarily indicate similarities in the underlying physical conditions.

To further illustrate this point, we mimicked the analysis on an observational set of star-forming complexes, treating each snapshot as an observation of a distinct cloud. We did not apply the biases of Sect. 5, to evince conclusions that do not depend on the assumptions made in that section. For each snapshot, we extracted the value of $\epsilon_{\text{ff}}$ corresponding to a given column density contour. We chose five different column densities as reference points: at 200, 500, 1000, 2000, and $5000 M_{\odot} \text{pc}^{-2}$. Consequently, we obtained four datasets (one for each feedback prescription) and computed the temporal average and dispersion of $\epsilon_{\text{ff}}$ for each simulation. The question is whether we can distinguish the different physical conditions of the clouds when we lack information about their temporal evolution and the dependence on $\Sigma_{\text{avg}}$.

The results, presented in Fig. 11, reveal that, except for NF, the averages of the different simulations are indistinguishable because of their significant temporal variations. If these were actual observations, the impact of stellar feedback would always be secondary compared to the variation due to the cloud evolution. This does not mean that feedback has no valuable impact on the SFR, since, looking at the full picture (Fig. 5), $\epsilon_{\text{ff}}$ changes substantially with $\Sigma_{\text{avg}}$ when the feedback prescription is changed. Therefore, the cause is that $\epsilon_{\text{ff}}$ is highly varying over
time and that the different runs achieve the same $\epsilon_{ff}$ during their evolution, but at different times. Only when looking at $\epsilon_{ff}(\Sigma_{\text{avg}})$ a clear dependence on the feedback mechanism becomes evident.

In addition, we notice that the differences just mentioned are mostly visible above the range of densities typically studied. Studies such as Lombardi et al. (2015), Pokhrel et al. (2021), or Bieging & Kong (2022) reach maximum column densities of the order of $1000 \Sigma_{\text{avg}}$, while the high resolution of these simulations allows us to probe an order of magnitude more. This density range is mainly governed by the cloud's structure, and observations may not be able to catch variations caused by feedback mechanisms.

Therefore, what is most interesting for observational datasets we can get with the present instrumentation are the differences in the slopes, which are detectable even at lower densities. From our study, we would expect clouds with strong radiation feedback to display the steeper slopes in the $\Sigma^* \Sigma_{\text{avg}}$ plot. It is extremely interesting to notice that in Pokhrel et al. (2021), clouds dominated by strong radiative feedback (Orion A, Sh2-140, Cep OB3, and Cygnus X) show on average a steeper slope than the other ones (Ophiucus, Orion B, Perseus, Aquila North, Aquila South, NGC 2264, AFGL 490, and Mon R2). In particular, the mean exponent of the KS relation found for the first clouds is 2.24 with a dispersion of 0.14, while the average for the second group is 1.87, with a dispersion of 0.24. Although the sample is limited, this supports our findings, and further studies that extend the sample could provide full confirmation of the trend.

## 7. Conclusions

We conducted an in-depth analysis of high-resolution magnetohydrodynamic simulations to investigate the impact of stellar feedback on two different star formation models: the original KS law and the Krumholz & McKee (2005) $\epsilon_{ff}$ model. Applying techniques commonly used in observational works, we reveal that substantial differences arise when stellar radiation and protostellar jets act to reshape the cloud structure. In particular, we find that the coupling between the area-column density relation and feedback mechanisms plays a major role in the observed star formation law, much more than the stellar distribution itself.

When only protostellar jets are active, they compress part of the gas such that the N-PDF becomes more weighted towards high densities. However, this does not translate directly into an increase in the SFR, instead resulting in a shallower KS law and a lower $\epsilon_{ff}$ at higher column densities. On the other hand, H II regions act to disperse the gas, leading to a steeper $\Sigma$ versus $\Sigma_{\text{avg}}$ relation, with an exponent that settles around a value of $q = 2$. This particular value also validates the study from Abreu-Vicente et al. (2015) conducted on Galactic molecular clouds. Even though the presence of ionising stars decreases the global SFR, the steepness of the $\Sigma(\Sigma_{\text{avg}})$ function is such that the observed star formation efficiency is higher at high column densities. The result is confirmed even for $\epsilon_{ff}$, although the weaker dependence on A mitigates the differences with the jet-only run.

We also applied prescription to emulate the effect of three different biases observations have to deal with. While assigning the same mass and age to each YSO detected does not significantly alter the observed SFR, the presence of a detection threshold for embedded YSOs and the limited resolution of column density maps can modify the shape of the relation at the highest densities. In particular, we show that low-resolution column density maps can affect the $\Sigma-\Sigma_{\text{avg}}$ relation when observing a highly structured cloud. In addition, we report a significant

impact due to the presence of a detection threshold, which leads to both a steep drop in the observed $\epsilon_{ff}(\Sigma_{\text{avg}})$ relation at high densities and an underestimation of the global SFR.

When applied altogether, the effect did not couple efficiently, so the resulting relation was close to that obtained by simply summing the biases. The only exception was NF, where the mass distribution of sinks changed towards the centre of the cloud. This caused an overestimation of the SFR due to the intertwined action of the YSO counting and the lower completeness of the detected embedded sources. Although this could be attributed to numerical effects, it could play a role when observing regions where the mass segregation mechanism is ongoing (either dynamical or due to the initial cluster formation). When all the biases are applied, we recover a good agreement with the observational results from Pokhrel et al. (2021), which suggests that the physics included in our simulation is sufficient to describe the process of star formation at cloud scales.

Finally, having access to the temporal evolution of our cloud, we could show that the exponent of the KS relation and the value of $\epsilon_{ff}$ can vary significantly with time. This means that in a sample of real clouds, where different evolutionary stages are present simultaneously, extracting the value of $\epsilon_{ff}$ at a single column density would hide the effect of feedback. Indeed, distinct departures between the models only become evident when considering the dependence on $\Sigma_{\text{avg}}$.

We will extend the sample of numerical simulations to cover a wider range of initial conditions and feedback strengths, to parameterise the impact of ionising stars as a function of cloud and stellar properties, and to perform a direct comparison with observational data.

**Acknowledgements.** We thank Mark R. Krumholz for the valuable insights provided. A.Z. thanks the support of the Institut Universitaire de France (IUF). We thank the referee for his relevant and focused insights, which have helped us enhancing the contextual relevance and clarity of the developments outlined in this work. This research has received funding from the European Research Council synergy grant ECOGAL (Grant: 855130).

## References

* Abreu-Vicente, J., Kainulainen, J., Stutz, A., Henning, T., & Beuther, H. 2015, A&A, 581, A74

* Allison, R. J., Goodwin, S. P., Parker, R. J., Portegies Zwart, S. F., & de Grijs, R. 2010, MNRAS, 407, 1098

* Antoniucci, S., Nisini, B., Giannini, T., & Lorenzetti, D. 2008, A&A, 479, 503

* Appel, S. M., Burkhart, B., Semenov, V. A., et al. 2023, ApJ, 954, 93

* Ascenso, J., Lada, C. J., Alves, J., Román-Zúñiga, C. G., & Lombardi, M. 2013, A&A, 549, A135

* Audit, E., & Hennebelle, P. 2005, A&A, 433, 1

* Barnes, A. T., Watkins, E. J., Meidt, S. E., et al. 2023, ApJ, 944, L22

* Bieging, J. H., & Kong, S. 2022, ApJ, 938, 145

* Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

* Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

* Bonnell, I. A. 2008, in ASP Conf. Ser., 390, Pathways Through an Eclectic Universe, eds. J. H. Knapen, T. J. Mahoney, & A. Vazdekis, 26

* Burkhart, B. 2018, ApJ, 863, 118

* Burkhart, B., Collins, D. C., & Lazarian, A. 2015, ApJ, 808, 48

* Burkhart, B., Stalpes, K., & Collins, D. C. 2017, ApJ, 834, L1

* Chen, H. H.-H., Burkhart, B., Goodman, A., & Collins, D. C. 2018, ApJ, 859, 162

* Chevance, M., Kruijssen, J. M. D., Krumholz, M. R., et al. 2022, MNRAS, 509, 272

* Clark, P. C., & Bonnell, I. A. 2005, MNRAS, 361, 2

* Dale, J. E., Ercolano, B., & Bonnell, I. A. 2013, MNRAS, 431, 1062

* Dobbs, C. L. 2015, MNRAS, 447, 3390

* Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11

* Elmegreen, B. G. 2011, in EAS Pub. Ser., 51, Star Formation in the Local Universe, eds. C. Charbonnel, & T. Montmerle, 45

* Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725

* Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

* Evans, Neal J., I., Heiderman, A., & Vutisalchavakul, N. 2014, ApJ, 782, 114
Federrath, C. 2013, MNRAS, 436, 3167

Federrath, C. 2015, MNRAS, 450, 4035

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43

Fiorellino, E., Tychoniec, Ł., Cruz-Sáenz de Miera, F., et al. 2023, ApJ, 944, 135

Frisch, U., & Bec, J. 2000, in New Trends in Turbulence, eds. M. Lesieur, A. Yaglom, & F. David (Springer EDP-Sciences), 341

Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

Grudić, M. Y., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 475, 3511

Gutermuth, R. A., Pipher, J. L., Megeath, S. T., et al. 2011, ApJ, 739, 84

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Harvey, P., Merín, B., Huard, T. L., et al. 2007, ApJ, 663, 1149

Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29

Hennebelle, P., Lebreuilly, U., Colman, T., et al. 2022, A&A, 668, A147

Hu, Z., Krumholz, M. R., Pokhrel, R., & Gutermuth, R. A. 2022, MNRAS, 511, 1431

Kennicutt, Robert C., J. 1989, ApJ, 344, 685

Kennicutt, Robert C., J., & De Los Reyes, M. A. C. 2021, ApJ, 908, 61

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Khullar, S., Krumholz, M. R., Federrath, C., & Cunningham, A. J. 2019, MNRAS, 488, 1407

Khullar, S., Federrath, C., Krumholz, M. R., & Matzner, C. D. 2021, MNRAS, 507, 4335

Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69

Kuhn, M. A., de Souza, R. S., Krone-Martins, A., et al. 2021, ApJS, 254, 33

Lada, C. J., Lombardi, M., Roman-Zuniga, C., Forbrich, J., & Alves, J. F. 2013, ApJ, 778, 133

Lada, C. J., Lewis, J. A., Lombardi, M., & Alves, J. 2017, A&A, 606, A100

Lee, Y.-N., & Hennebelle, P. 2018, A&A, 611, A89

Lee, E. J., Miville-Deschênes, M.-A., & Murray, N. W. 2016, ApJ, 833, 229

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2017, ApJ, 846, 71

Lombardi, M., Lada, C. J., & Alves, J. 2013, A&A, 559, A90

Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45

Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, A1

Ma, Y., Wang, H., Zhang, M., et al. 2022, ApJS, 262, 16

Mayya, Y. D., Alzate, J. A., Lomelí-Núez, L., et al. 2023, MNRAS, 521, 5492

Mazumdar, P., Wyrowski, F., Urquhart, J. S., et al. 2021, A&A, 656, A101

Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5

Megeath, S. T., Gutermuth, R. A., & Kounkel, M. A. 2022, PASP, 134, 042001

Menon, S. H., Federrath, C., & Kuiper, R. 2020, MNRAS, 493, 4643

Menon, S. H., Federrath, C., Klaassen, P., Kuiper, R., & Reiter, M. 2021, MNRAS, 500, 1721

Murray, N. 2011, ApJ, 729, 133

Murray, D., Goyal, S., & Chang, P. 2018, MNRAS, 475, 1023

Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41

Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40

Palmeirim, P., Zavagno, A., Elia, D., et al. 2017, A&A, 605, A35

Pokhrel, R., Gutermuth, R. A., Krumholz, M. R., et al. 2021, ApJ, 912, L19

Pomarès, M., Zavagno, A., Deharveng, L., et al. 2009, A&A, 494, 987

Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431

Porter, D. H., Pouquet, A., & Woodward, P. R. 1994, Phys. Fluids, 6, 2133

Raboud, D., & Mermilliod, J. C. 1998, A&A, 333, 897

Roccatagliata, V., Preibisch, T., Ratzka, T., & Gaczkowski, B. 2013, A&A, 554, A6

Salpeter, E. E. 1955, ApJ, 121, 161

Sandford, M. T., I., Whitaker, R. W., & Klein, R. I. 1982, ApJ, 260, 183

Schmidt, M. 1959, ApJ, 129, 243

Schmidt, M. 1963, ApJ, 137, 758

Schneider, N., Ossenkopf, V., Csengeri, T., et al. 2015, A&A, 575, A79

Seifried, D., Walch, S., Girichidis, P., et al. 2017, MNRAS, 472, 4797

Spilker, A., Kainulainen, J., & Orkisz, J. 2021, A&A, 653, A63

Stolte, A., Grebel, E. K., Brandner, W., & Figer, D. F. 2002, A&A, 394, 459

Suin, P., Shore, S. N., & Pavlík, V. 2022, A&A, 667, A69

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2023, ApJ, 945, L19

Tassis, K., Christie, D. A., Urban, A., et al. 2010, MNRAS, 408, 1089

Teyssier, R. 2002, A&A, 385, 337

Usero, A., Leroy, A. K., Walter, F., et al. 2015, AJ, 150, 115

Utomo, D., Sun, J., Leroy, A. K., et al. 2018, ApJ, 861, L18

Vacca, W. D., Garmany, C. D., & Shull, J. M. 1996, ApJ, 460, 914

Verliat, A., Hennebelle, P., González, M., Lee, Y.-N., & Geen, S. 2022, A&A, 663, A6

Wall, J. E., Mac Low, M.-M., McMillan, S. L. W., et al. 2020, ApJ, 904, 192

Willis, S., Guzman, A., Marengo, M., et al. 2015, ApJ, 809, 87

Zari, E., Lombardi, M., Alves, J., Lada, C. J., & Bouy, H. 2016, A&A, 587, A106

Zhang, S., Zavagno, A., Yuan, J., et al. 2020, A&A, 637, A40

Zhang, S., Zavagno, A., López-Sepulcre, A., et al. 2021, A&A, 646, A25
**Appendix A: Effect of projection axis**

In our analysis we mainly show the results from the projection
along the z-axis. However, the other axes behave similarly and no
notable difference has to be reported. Figure A.1 shows how the

KS relation changes when averaging between the three lines of
sight. The shaded areas represent the 1 σ dispersion of the three
curves. It is visible that the variability decreases with time, and
the results remain coherent. This could also be guessed looking
at Fig. 6, where the dispersion of the fitted slope is small.

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Same as Fig. 5, but using the average and 1 σ dispersion of ϵ$_{ff}$ obtained from the three lines of sight. Data are masked when at least one of the three contours at the same density threshold does not satisfy the 500 pixel criterion.</figcaption></figure>
## Appendix B: 2D versus 3D star formation relation

The non-trivial correlation between local quantities and their projected value constitutes a major challenge to the interpretation of astrophysical observations. With simulations, one can accurately preserve the evolution of physical variables for each cell, accessing the full 3D information. Therefore, it is of great interest to understand the extent to which projection effects impact the results. Hence, we extended the same cumulative method described in Sect. 3 to perform the 3D analysis. We divided the cloud into 500 volume density contours, evenly spaced in log between $10^3$ and $10^7$ cm$^{-3}$. Converting the 500 pixel criterion used for 2D contours, we require the isosurfaces to contain at least a volume greater than $500^{3/2} \cdot (7.3 \times 10^{-3} \text{pc})^3 \approx 4 \times 10^{-3} \text{pc}^3$. We plotted the outcome in Fig. B.1, where now we scaled the y-axis (the volume SFR density $\dot{\rho}_{\star}$) by $\rho_{\text{gas}}^{1.5}$, where $\rho_{\text{gas}}$ is the average volume density of the contour. The choice of the exponent is justified by Fig. 5. Since Eq. (1) displays a constant behaviour with $\Sigma_{\text{avg}}$ in the radiative simulations, we expect a similar feature when dividing by $\rho_{\text{gas}}^{1.5} \propto \rho_{\text{gas}}/t_{\text{ff}}$.

The first thing we could estimate from the plots is the limit up to which data can be considered reliable. Indeed, we expect the efficiency to increase sharply as we approach the density threshold for sink formation (Khullar et al. 2021), and we can see a steep rise towards the very high-density edge of the graph. Keeping in mind that caution is always necessary when comparing 3D and 2D datasets, the distance between the knee of the rise and our reliability limit (beginning of the dashed line) gave us confidence that our constraint provides sufficient protection against resolution effects.

Comparing these graphs with their respective counterparts in Fig. 4, we can confirm that the features observed in the projection are supported by the 3D analysis. The plot at 2.1 Myr shows that its 3D version closely resembles the 2D version. At $t = 2.7$ Myr, runs with jets exhibit similar trends. At the same time, the radiation-only simulation presents a comparable $\dot{\rho}_{\star}$ to the run with both jets and radiation at low densities, but an higher efficiency at higher densities. Graphs at 3.0 and 3.5 Myr, instead, require more attention. Indeed, the snapshot at 3.0 Myr confirms the findings of the previous section, with the relative behaviours that remain preserved, even if the differences appear to be smaller. In particular, one may notice that the departure of HIIR+PSJ from PSJ happens ‘later’. However, comparing points between the two graphs is dangerous because, without a characteristic scale, there is no exact conversion from column to volume density.

At 3.5, Myr, the slope displayed by simulations with radiation flattened out even more. However, there are some notable differences. In the HIIR run, the star formation event has nearly ceased, and the SFR density curve lies below the others. On the other hand, the HIIR+PSJ run achieves a $\dot{\rho}_{\star}$ comparable to PSJ at high column densities $n$, although this happens just below the reliability limit. Nevertheless, this does not affect the conclusions drawn so far. Indeed, at other times in the simulations, the observed trends are consistent with the previous discussions.

<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Evolution of $\dot{\rho}_{\star}$ as a function of volume number density. The dashed lines mark the 3D contours with volumes smaller than $4 \times 10^{-3} \text{pc}^3$. Units of the y-axes are $M_{\odot}^{-0.5} \text{pc}^{1.5} \text{Myr}^{-1}$. The steep rise to the right is caused by the density threshold for sink formation at $n = 10^7 \text{cm}^{-3}$.</figcaption></figure>