Star-formation-rate estimates from water emission

(Abridged) The star-formation rate (SFR) quantitatively describes the star-formation process in galaxies. Current ways to calibrate this rate do not usually employ observational methods accounting for the low-mass end of stellar populations as their signatures are too weak. Accessing the bulk of protostellar activity within galactic star-forming regions can be achieved by tracing signposts of ongoing star formation. One such signpost is molecular outflows, which are bright in molecular emission. We propose to utilize the protostellar outflow emission as a tracer of the SFR. In this work, we introduce a novel version of the galaxy-in-a-box model, which can be used to relate molecular emission from star formation in galaxies with the SFR. We measured the predicted para-H2O emission at 988 GHz and corresponding SFRs for galaxies with LFIR = $10^8$ - $10^{11}$ L$_\odot$ in a distance-independent manner, and compared them with expectations from observations. We evaluated the derived results by varying the star formation efficiency, the free-fall time scaling factor, and the initial mass function. For the chosen H2O transition, relying on the current Galactic observations and star formation properties, we are underestimating the total galactic emission, while overestimating the SFRs, particularly for more starburst-like configurations. The current version of the galaxy-in-a-box model accounts for a limited number of processes and configurations, that is, it focuses on ongoing star formation in massive young clusters in a spiral galaxy. Therefore, the inferred results, which underestimate the emission and overestimate the SFR, are not surprising: known sources of emission are not included in the model. To improve the results, the next version of the model needs to include a more detailed treatment of the entire galactic ecosystem and other processes that would contribute to the emission.


Introduction
Star formation lies at the very center of the baryon cycle and plays a pivotal role in shaping galactic ecosystems. There are different measures of this process, all of which help to understand and characterize its behavior through cosmic history. One example is the star-formation rate (SFR), as it provides a quantitative description of the star-forming properties of a given object by relating the total mass of stars formed in a give time unit, that is, M * /∆t. The SFR is used to establish the cosmic star formation history (e.g., Lilly et al. 2013;Madau & Dickinson 2014), which in turn is used to understand and quantify the evolution of galaxies.
The key epoch of cosmic star formation history, which is when star formation peaked (known as cosmic noon), marks a critical stage during the evolution of today's galaxy population (e.g., Shapley 2011;Madau & Dickinson 2014;Förster Schreiber & Wuyts 2020). Cosmic-noon galaxies, lying at redshifts of 2-3, exhibit extremely different SFRs from those observed in the local Universe, reaching >1000 M yr −1 , while the Milky Way is forming stars at a rate of ∼1 M yr −1 (e.g., Kennicutt & Evans 2012).
There are various ways of deriving the SFRs in galaxies from nebular line, UV, infrared, radio, and X-ray emission (Madau & Dickinson 2014). These methods all assume that there is a scaling between the luminosity in a given band and the SFR. However, the observed emission is usually dominated by high-mass stars, which easily outshine low-mass stars due to their energetic output, and so an initial mass function is applied to correct for low-mass stars, which is where most of the mass resides. In the local Universe, the SFR is readily traced and calibrated with Hα, Hβ, [O ii], and [O iii] emission (e.g., Kennicutt 1998;Tresse et al. 2002;Kewley et al. 2004;Salim et al. 2007;Villa-Vélez et al. 2021). However, in the past 20 years, advances in astrochemistry have provided additional ways to trace star formation, even in its most embedded stages, and allow us to trace the episodes of current star formation in galaxies (e.g., Herbst & van Dishoeck 2009;Jørgensen et al. 2020).
Molecular emission from protostars is not yet commonly used as a SFR tracer. Nevertheless, this emission has the potential to trace even low-mass populations directly. At the earliest stages, the forming star itself is deeply embedded in gas and dust and is thus completely obscured. Therefore, the key is to trace signposts of these early stages that are not obscured. One of these signposts is outflows, which are launched from protostars in their main accretion phase when the interaction between the infalling envelope, winds, and jets launched from the protostar are particularly strong (Bally 2016). These outflows are launched from close to the protostar, but quickly punch their way through to the surrounding molecular cloud, where they are not obscured (Bachiller et al. 1990). In our Galaxy, one of the best tracers of this protostellar component is water (van Dishoeck et al. 2021), which is predominantly locked up as ice on dust grains, but is released from the grain mantles into the gas phase, causing a jump in the abundance of many orders of magnitude. At the same time, the physical conditions are conducive to water being readily excited into rotational states (e.g., Suutarinen et al. 2014).
Water emission is also observed toward high-redshift galaxies (e.g., Yang et al. 2013Yang et al. , 2016Jarugula et al. 2019), where it too has been calibrated to serve as an SFR tracer (Jarugula et al. 2019). However, at high redshift, water is thought to trace dusty molecular clouds illuminated by either massive stars or a central galactic nucleus, and therefore the excitation is assumed to be via far-infrared (FIR) pumping (e.g., González-Alfonso et al. 2008. However, toward the Galactic sources, which were extensively observed with the Herschel Space Observatory (e.g., the Water In Star-forming regions with Herschel survey (WISH); van Dishoeck et al. 2011Dishoeck et al. , 2021, water emission is almost uniquely associated with outflows, where its excitation is collisionally dominated, and other processes, such as FIR pumping, have a negligible contribution to the excitation Goicoechea et al. 2015).
With the goal of tracing active star formation in galaxies with molecular emission from protostars,  created a galactic model, the so-called galaxy-in-box model, simulating emission from star-forming regions. Using up-to-date knowledge of Galactic star formation and state-ofthe-art astrochemical observations of Galactic protostars, the galaxy-in-a-box model simulates emission from young clusters in a chosen galaxy, and at the same time provides insight into the statistics of the star formation process. The default molecular emission is that from water at 988 GHz (J Ka,Kc = 2 02 − 1 11 ), which is readily observed even at high redshifts where its emission is thought to be dominated by the FIR-pumping-dominated regions, as outlined above.
In this work, we present an extension to the galaxy-in-a-box model, which allows us to derive SFRs from simulated galaxies and their individual star-forming clusters, and to put constraints on local and global SFRs. We focus on water emission at 988 GHz, and simulate emission for galaxies with L FIR = 10 8 −10 11 L for varying star-formation parameters. This paper is organized as follows. Section 2 describes all of the changes introduced to the galaxy-in-a-box model. Subsequently, in Sect. 3 we present the results of this study and test them against observations and the literature; we then discuss these comparisons in Sect. 4. Finally, we present our conclusions in Sect. 5.

Methods
In this study, we explore the relation between the SFR, water luminosity (L H 2 O ), and far-infrared luminosity (L FIR ) using the galaxy-in-a-box model , for an overview of the model see Appendix A). This is a novel, stateof-the-art astrophysical modeling tool that simulates emission from young clusters in a galaxy and provides detailed insights into the constituents of the star-formation process and derived parameters. The model relies on relatively few input parameters, giving the user great flexibility to define global and local galactic parameters.
For deriving the SFR and relating L FIR to the virial mass of galaxies, we implemented a number of upgrades to the model, which we describe in Sect. 2.1. In Sect. 2.2 we describe the choice of parameters for the simulated galaxies.

Changes to the galaxy-in-a-box model
For the purposes of this study, we introduced the SFR as an input and output parameter in the galaxy-in-a-box model. We only used the output SFRs. However, in the following, we describe the full extent of the new SFR feature. The SFR tells us how much material is turned into stars per unit of time. With that in mind, we defined the SFR for a cloud in a galaxy as where N * is the number of formed protostars, M * is the average protostellar mass, t cloud is the age of the cloud, τ sc ff is the unitless free-fall scaling factor , and t ff is the free-fall time of the cloud. In the case of the galaxy-ina-box model, the age is randomized, that is, it randomly scales the ages such that they range from newly formed to completely collapsed. The global SFR of the entire galaxy is therefore the sum of the individual rates for each cloud.
In the model, we assume that each cloud goes on to form one cluster; in nature, clouds may go on to form several generations of clusters, but for the purposes of this study, where we consider global star formation, this is not relevant. With this implementation of the SFR, we introduce the possibility to also constrain the SFR at the cloud or cluster level. The cluster module can now be run with a fixed SFR, where the age of the cluster is adjusted through the free-fall time scaling factor, which can be easily derived from Eq. (1): In this equation, t ff is already randomized (t ff,random ) to avoid poor SFR adjustment due to age assignment that takes place later in the model. However, Eq. (2) is not used in this study. On a global scale, that is, when introducing constraints on the total galactic SFR, the new version of the galaxy-in-a-box model monitors the total SFR of the given galaxy and computations stop when the specified SFR is reached. The allowed deviation from the specified SFR is ±10%. There may be situations where the galaxy-in-a-box model will not converge: these situations are unphysical, and an example would be a very low-mass galaxy with a very high SFR. There needs to be enough gas that can be turned into stars at the desired rate.
One of the changes to the galaxy-in-a-box model that was essential for this study was to set the number of clusters as limited by the total molecular gas reservoir, rather than setting it as a fixed number in the input file. We infer the number of molecular clouds from which star forming clusters form by putting an upper limit on the total mass of the molecular clouds, which are randomly generated using the molecular cloud mass distribution. When the limit is reached, the clouds are no longer passed to the cluster part of the calculations (see Fig. 1 of . This way we ensure that the mass of clouds does not exceed the available molecular reservoir. Lastly, the mass properties of the galaxy can now be set by defining the L FIR of the galaxy.  Good (1989). The solid straight line represents the best-fit power law to the data points, the darker shaded region corresponds to the 95% confidence region of the correlation, and the lighter shaded region represents the region that contains 95% of the measurements.
through the observed M vir − L IR relation (see Fig. 1). The viral mass of the galaxy can be expressed as: where L IR is the total far-infrared luminosity of the cloud. With Eq.
(3), we can simulate L H 2 O for galaxies with different L IR , including typical galaxy types observed with H 2 O emission, that is, subluminous infrared galaxies (subLIRGs; L IR < 10 11 L ), LIRGs (10 11 L ≤ L IR < 10 12 L ), ultraluminous infrared galaxies (ULIRGs; 10 12 L ≤ L IR < 10 13 L ), and hyperluminous infrared galaxies (HyLIRGs; L IR ≥ 10 13 L ). In this study, we are interested in relative values for derived luminosities and SFRs, and therefore we make an assumption that L FIR is a proxy for L IR , and use these luminosities interchangeably.

Considered parameters
The goal of this study is to explore the SFRs derived with the galaxy-in-a-box model and how they relate to derived luminosities. To achieve this goal, we decided to use the template galaxy from  with emission from the para-H 2 O 2 02 − 1 11 line at 987.927 GHz, and tweak the star formation efficiency, the free-fall-time scaling factor, and the initial mass function. The exact ranges of the probed parameters are described in Table 1. For the galactic masses, or in this case luminosities, we decided to probe galaxies with L FIR = 10 8 −10 11 L , where for the range 10 8 −10 10 L we continued with an increment corresponding to the given order of magnitude (i.e., 10 8 , 2 × 10 8 , 3 × 10 8 , etc.), and we stopped at 10 11 L . We made this choice because we wanted to probe the chosen regime in a relatively uniform way. Moreover, the lower limit was dictated by low galactic mass (10 8 L corresponds to ∼10 7 M ), while the upper one was dictated by the limitations of the computational power. As shown below, the inferred SFRs can readily be extrapolated to even higher luminosities. Table 1. Parameters considered in this study.
Red filling corresponds to parameter space used for the high-z correlation test, while x refers to those considered together with the standard correlation. For the latter, the only omitted galactic type was that with L FIR = 10 10 L for all combinations involving ε SF = 30%. In the IMF parameters, s refers to standard, t-h to top-heavy, and b-h to bottomheavy.
In the model, we use the relation between mass and H 2 O line luminosity obtained only from Galactic sources to estimate the amount of emission generated by protostars. As a sanity check, we can include the high-z observations in this correlation, as shown in Fig. 2. Including the high-z measurements shifts the correlation slightly, such that low-mass protostars are assigned less emission, and vice versa for high-mass protostars. Therefore, if these high-z sources are included, this has implications for the assumed IMF. In order to obtain luminosity distances for high-z objects, we used a Planck 2018 flat ΛCDM cosmology with H 0 = 67.7 km s −1 Mpc −1 and Ω M = 0.310, as implemented in the Astropy package (Astropy Collaboration 2018).

Results
By extracting SFRs together with L H 2 O , while at the same time defining galaxies according to their luminosity rather than their mass directly, we are able to confront expectations based on the literature about the star formation process as seen in the Milky Way while simultaneously testing the galaxy-in-a-box model. Therefore, in this proof-of-concept study, we ran a number of simulations spanning a range of parameters representing different galactic and star formation properties (see Table 1).
As mentioned in Sect. 2.1, we used two different mass-line luminosity correlations: in the first, we only use the Galactic data points, and in the second we include the high-z data points. We excluded certain parameters from the high-z test, because they were either computationally heavy or unnecessary for testing the impact of the high-z extrapolation (for further discussion, see Sect. 4). For galaxies with L FIR = 10 8 −9 × 10 8 L , we ran 40 simulations for each setup, while for other luminosity ranges we ran 20 simulations per setup. The increased number of simulations for this specific galactic type was dictated by higher SFR variations, as the molecular reservoir is relatively low, which is reflected in larger variations in the number of formed stars. We also excluded calculations for galaxies with L FIR = 10 10 −9×10 10 L , which would have ε SF = 30%, because they were the most computationally heavy, and including them would not affect any conclusions of this study. In total, we ran 15 200 simulations, including 12 240 main runs and 2960 runs for the high-z test. Uncertainties for the simulation results are calculated as a standard deviation from the mean value, which is derived for all runs with the same set of parameters. The best fits were obtained using linear regression while accounting for the spread in the y-direction. If the spread is not shown, this means that the size is smaller than the marker or the line size. When recalculating fluxes to luminosities, we naturally account for the propagation of uncertainties.
We describe the literature sample chosen for this study in Sect. 3.1. We then present the results through derived Sect. 3.4), and L H 2 O − SFR (Sect. 3.5) relations, which we compare to those provided in the literature.

Literature sample
As a default source of Galactic observations, we use data from the Water Emission Database   . This sample includes both nearby subLIRGs, LIRGs, and quasars, as well as high-z quasars, ULIRGS, and HyLIRGs, with the farthest one being the HyLIRG, namely HFLS3, at z = 6.337 (D L = 62834.75 Mpc; for more details see Riechers et al. 2013). A detailed description of the sample and exact values used in this study can be found in Kristensen et al. (2022).

Total stellar mass versus SFR
We evaluated the derived SFRs by exploring their relation with the total stellar mass of the corresponding galaxies. From Fig. 3, we see that we are overestimating the SFRs when looking at functions derived by for example Salmon et al. (2015) for the main sequence galaxies and Rinaldi et al. (2022) for the starbursts.
With the chosen set of properties, galaxies with M * < 10 6.5 M seem to lie close to the main sequence estimates from Salmon et al. (2015), at least in their lower limits. However, going to cases where the combination of considered parameters resulted in an increase in SFRs, especially galaxies with M * > 10 6.5 M , we start overestimating SFRs by at least one order of magnitude when compared to the literature (Rinaldi et al. 2022).
We also observe two distinct populations that appear to be dictated by the value of the free-fall-time scaling factor. For τ sc ff = 1, we let the efficiency of the free-fall time depend only on the density of the progenitor molecular cloud, while by introducing τ sc ff = 5 we prolong the time required to form most of the stellar population, resulting in a more diverse range of protostellar ages. From Fig. 3, we see how a decrease in the free-fall-time scaling factor influences the derived SFR. Considering the relatively low efficiency of the star formation process, the lower-SFR branch is likely to be more consistent with the nature of the star formation process. We discuss this topic further in Sect. 4.4.

L FIR − L H 2 O correlation
To compare the predicted fluxes with observations, we first converted them to luminosities using the following expression: where I is the total intensity in Jy km s −1 , λ 0 the wavelength in microns (303.4557 µm for the para-H 2 O 2 02 − 1 11 line), and D L the luminosity distance of the source in megaparsecs. By converting fluxes, we can quantitatively compare our results with observations, as they are no longer distance dependent. Using linear regression, we derived best-fit lines to the following expression: log 10 L H 2 O /L = a × log 10 (L FIR /L ) + b.
Table 2 provides all of the derived slopes and intercepts. In the following, we focus on the two setups exhibiting the highest and lowest water emission. These are the models with ε SF = 30%, IMF = top-heavy, τ sc ff = 1, and ε SF = 1%, IMF = bottom-heavy, τ sc ff = 5, respectively. For the least emitting case, we derive a = 0.809 ± 0.003 and b = −7.269 ± 0.029, while for the most emitting case we derive a = 0.809 ± 0.001 and b = −5.135 ± 0.012. In both cases, R 2 = 99.9%. For all of the simulations, the slope stays roughly constant with a ≈ 0.81, and therefore the span of luminosities is described by the intercept falling in the range of −7.269 to −5.135. We can derive the general relation for waterline luminosity depending on the intercept value: From Fig. 4, we see that we deviate from extragalactic observations by between a factor of a few and about two orders of magnitude. We observe that the expectations built on the extragalactic sample taken from Jarugula et al.   Fig. 3. SFR as a function of stellar mass of each galaxy. Full color markers represent results, where the free-fall-time scaling factor was set to 1, while markers in the same but lighter colors correspond to τ sc ff of 5. Circles represent setups with the standard IMF (Chabrier 2003), while triangles pointing upwards and downwards represent setups with its top-heavy and bottom-heavy versions, respectively. Different colors of the markers refer to different star formation efficiencies, where green, orange, and red mean an ε SF of 1%, 10%, and 30%, respectively. Solid lines represent best-fit lines from Rinaldi et al. (2022) to their starburst (SB) population, while dashed lines represent best fits to the main sequence galaxies from Salmon et al. (2015). in Sect. 4.1) -are especially far from our expectations for the brightest high-z galaxies. We discuss this further in Sect. 4.5. Also, in Sect. 4.2, we explore the possible impact of the inclusion of high-z starbursts on the correlation between the envelope mass and intensity (M env −I relation) -which is the basis of emission assignment in the galaxy-in-a-box model -and whether it could explain the observed differences.

L FIR − SFR correlation
To further evaluate derived SFRs, we explored their relation with corresponding far-infrared luminosities (Fig. 5). We clearly see that the derived SFRs create different populations depending on the star formation efficiency and the free-fall-time scaling factor. Again, we are clearly overestimating the SFRs. However, relations in the literature, for example those of Kennicutt & Evans (2012) and Casey et al. (2014), fall into our lower prediction regime, meaning that at least for the star forming galaxies with lower star formation activity (with respect to the standard setup in the galaxy-in-a-box model), we are roughly recovering the expected star formation process.
The span of the SFRs derived in this study depends strongly on the efficiency of the process. The discrepancy between the literature values and our simulations can be as high as two orders of magnitude. We focused on and derived relations analogous to Eq. (5) for the setups with the lowest and highest emission, as well as the standard model setup from the galaxy-in-a-box model. We provide all of the derived relations in Table 2. Here, we do not derive almost identical slopes, as we did for L FIR − L H 2 O . For the most extreme cases of the L FIR -SFR relation, we derive slopes of 0.94 ± 0.04 and 0.90 ± 0.03, which agree within the uncertainties, while the derived intercepts (here, the intercept refers to the term b in Eq. (5), which is further used as showed in Eq. (6)) are equal to −8.50 ± 0.35 and −5.75 ± 0.33, respectively. We further discuss the apparent excess in SFR in Sect. 4.3.

L H 2 O − SFR correlation
The last explored dependence was that of L H 2 O and the corresponding SFRs. We see from Fig. 6 that all of the derived SFRs fall into the same population, which is expected considering the fact that the greater the luminosity, the more actively starforming and massive the corresponding galaxy. By fitting all of the derived points to Eq. (5), we get a slope of 1.11 ± 0.01 and an intercept of −0.083 ± 0.018, indicating a near-proportionality between the SFR and L H 2 O .

Galactic type
Model results log 10 Y = a · log 10 X + b Notes. Results from running simulations with all considered parameters combinations (for more details see Table 1). In the IMF parameters, s refers to standard, t-h to top-heavy, and b-h to bottom-heavy. For details on the τ sc ff and its relation with the free-fall time efficiency we refer the reader to , where extensively discussed the impact of this factor. The four setups in the bottom of the table divided by the horizontal line refer to our high-z test (see Sect. 4.2). Since the test had impact only on the water emission, we did not provide the SFR -L FIR relations as these are the same as for the corresponding Galactic setups. The correlation between the L H 2 O -L FIR had consistent R 2 ≈ 0.999, hence we do not provide the R 2 values for this relation in the table.
However, Fig. 6 suggests that we are systematically overestimating SFRs by approximately four orders of magnitude with respect to the findings of Jarugula et al. (2019), where SFR M yr −1 = 7.35 +5.74 −3.22 ×10 −6 L H 2 O (L ). If extrapolating their relation to Galactic star-forming regions, we would underestimate SFRs by orders of magnitude . We discuss this discrepancy in Sect. 4.5.

Discussion
In the following, we discuss derived SFRs and water luminosities. We also evaluate how the star-formation parameters considered here could affect the results and compare our results with the literature. Moreover, we discuss what other physical processes not considered in this study could impact the derived values and explore other possible influences.

Insights from L H 2 O /L FIR ratios
The ratio of L H 2 O and corresponding L FIR could be used to understand the source of the observed water emission (this is shown in Fig. 7). This in turn can help us to understand whether or not water behaves differently in different galactic regions and galactic types. With this in mind, we calculated the ratios derived from the galaxy-in-a-box model and compared them with our Galactic and extragalactic samples.
The derived values 10 −8 < L H 2 O /L FIR < 10 −6 fall below those from all objects considered in the extragalactic sample, but coincide with the Galactic sample at its high-mass/highluminosity end (see Fig. 4). We know from Galactic observations (e.g., van Dishoeck et al. 2021) that water emission from young stellar objects predominantly comes from the shocked material in outflows. Therefore, a natural assumption would be that the Galactic sample is consistent in terms of the calculated ratios. Instead, what we see is that low-to intermediate-mass protostars exhibit roughly the same ratios as the extragalactic sample, and we see a clear drop for the most luminous end of the Galactic objects.
Available water observations of Galactic high-mass young stellar objects are limited because of both their number and sensitivity. One of the most detailed studies was conducted with a survey towards the Cygnus X star-forming region (PI: Bontemps; San José-García 2015). Cygnus-X is one of the nearest massive star-forming complexes (D ∼ 1.3-1.4 kpc, e.g., Rygl et al. 2012). However, even these observations do not recover the total emission that would come from a highmass-star-forming complex because of the spatial resolution and sensitivity limitations of the HIFI instrument on the Herschel Space Observatory. This latter survey, one of the most complete, only consists of single-pointing observations. The gray-shaded area between these two lines refers to the probed parameter space, and all possible outcomes considered in this study would fall in that regime. Solid blue and red lines refer to the results derived for setups with the top-heavy IMF form for the Galactic and extragalactic M env − I relations, respectively. Dotted lines show the results for these two correlations, when the standard IMF is applied. In both, i.e., the standard and the top-heavy cases, the free-fall-time scaling factor is set to 1. Circles refer to observational samples (for more details we refer the reader to Sect. 3.1), while the purple line represents the expected relation from Jarugula et al. (2019). Therefore, new instruments are needed to fully estimate the amount of H 2 O emission coming from a forming Galactic cluster.
To take another approach, we estimate the amount of H 2 O emission from the nearby W3 high-mass-star-forming region. Its distance is 2 kpc and its age is 2 Myr (Bik et al. 2012). We used a mass of 4×10 5 M for the entire cluster (Rivera-Ingraham et al. 2013), corresponding to a total luminosity of 2×10 6 L using Eq. (3). To estimate the missing emission from all protostars, we ran a model for just one cluster instead of an entire galaxy. The cluster model predicts a total line intensity of 120 K km s −1 , which may be compared to the observed value of the highmass protostar W3-IRS5 of 21.9 K km s −1 (van der Tak et al. 2013), which has a luminosity of 10 5 L , or 5% of that of the cluster. The simulated value is highly sensitive to the adopted age of the cluster, for example, 1 Myr would result in a predicted intensity of 250 K km s −1 . This implies that for an individual cluster, we need to know the age accurately to within 10%, which is not currently possible. It is reasonably possible that the amount of water emission we are missing is between a factor of 6 and 12. Without being able to map the entire cluster in water emission, we will not know exactly how much.

High-z test
Knowing that the relation between water emission and L FIR spans over many orders of magnitude starting from the low-mass protostars to high-z HyLIRGs, we probed the influence of the extragalactic observations on the M env − I relation, and explore how this extrapolated form of the formula impacts the derived intensities.
In Fig. 2, we see that by including the extragalactic observations, we effectively lower the contribution from the low-mass end of the correlation and we see that it will only positively impact the high-mass protostars. On the other hand, the purely Galactic correlation lowers the emission from the high-mass protostars. Therefore, considering that we are underestimating water emission, we focused only on the standard and top-heavy IMF forms. We did this because the standard IMF is already dominated by the low-mass end of the distribution, and we also know from  that the emission derived for the bottom-heavy IMF is practically indistinguishable from the standard one. At the same time, the top-heavy IMF would increase the emission even for the normal form of the correlation, and the inclusion of the extragalactic sources increases the slope by ∼10% (see Fig. 2). The results of the test indicate that inclusion of the extragalactic sources results in lowered emission, on average, and that the difference with the results with the purely Galactic correlation starts to diminish for higher galactic masses and higher star-formation efficiencies. This effect is not surprising as the star-formation process is dominated in both total mass and number by low-mass protostars, while in terms of total bolometric luminosity the high-mass stars completely dominate the picture (e.g., Kroupa 2002). Therefore, the inclusion of the extragalactic sources, which lowers the emission from the low-mass protostars, naturally lowers the water emission derived from the simulated galaxies, as this is the main star-forming component if we consider Milky Way-like star formation. However, for the high-z starbursts with high star-formation efficiencies and seemingly top-or even extremely top-heavy IMFs, this extrapolation could make a difference, when simulating star formation and its emission. Nevertheless, we do not investigate this further, as this is beyond the scope of this paper.

SFR estimates
From the results derived in this study, we are consistently overestimating SFRs for given galactic types. However, when considering the assumptions behind the model, and the fact that in the current version of the model, current SFRs are simulated without correcting for star formation histories or existing populations, the overestimation is no longer prominent. The galaxy-in-a-box model was created as a tool for simulating emission from active and current star formation in galaxies. Therefore, even though the model accounts for dynamical differentiation of (proto)stellar ages, the model does not account for already existing, older stellar populations that normally would contribute to observations from which the rates are calculated. Moreover, as seen in Figs. 3 and 5, the results lie close to the literature estimates, if we assume low star formation activity. A calibration of the SFRs of galaxies depends on their current star formation activity. If the bulk of galaxies are observed during a period of low star formation, we would naturally fall on the lower SFR side. Also, there are many factors influencing star formation activity in galaxies that are not taken into account in the current version of the galaxy-in-a-box model.
Another important aspect is that when calibrating SFRs from L FIR , one has to make assumptions about parameters such as the IMF and star formation history, which are sources of additional uncertainty in the final estimation of the SFR. Moreover, L FIR is likely to underestimate the SFR in young clusters (Gutermuth et al. 2011) by up to an order of magnitude, and these are the main objects of interest in this study. If this is the case, our SFR estimates are roughly consistent with expectations.
Lastly the galaxy-in-a-box model accounts for all stellar products, from brown dwarfs to high-mass stars. Therefore, it is not subject to observational limitations and the apparent overestimation could be an effect of accounting for all objects, including those that are normally unobservable, as illustrated in the W3 example above. The scenario we are considering slightly more closely resembles the high-z situation, where galaxies are filled with active star-forming regions and are described as 'full of Orions' (Rybak et al. 2020). In this case, having relatively young star-forming regions, we trace only active and current star formation without accounting for higher differentiation of ages and stellar populations.

Impact of the star-formation parameters
In this study, we explored simulations for different galaxy types, and as such explored a broad parameter space (see Sect. 2.2 and Table 1). As in the first galaxy-in-a-box study , we observe no strong effect of the IMF, even though we included nearby subLIRGs, LIRGs, and quasars, as well as high-z quasars, ULIRGS, and HyLIRGs in the correlation that is used to assign molecular emission to protostars. This is expected as the extrapolation to the high-z regime changes the slope of the correlation only by ∼10%.
We observe a strong impact of the star-formation efficiency and the free-fall-time scaling factor, both for the derived emission and SFRs. This is of no surprise as both parameters impact the stellar population of each cluster. The free-fall-time scaling factor will effectively lower the ages of the clouds and thus increase the emission, while the star-formation efficiency regulates how much of the molecular reservoir will be turned into stars, hence increasing the number of stars.
One of the new input parameters in the galaxy-in-a-box model is the mass of the galaxy, as derived from Eq. (3). Clearly, the more massive the galaxy, the more emission we derive from the model. However, this parameter has its own uncertainty, which would be especially important when considering the predicted water emission. The relation between the mass and luminosity was also derived for young stellar objects by Pitts et al. (2022), where: log (M env /M ) = 0.30 +0.07 −0.06 + 0.79 +0.01 −0.02 log (L bol /L ) .
Although this expression was inferred for individual protostellar envelopes, it clearly agrees with Eq. (3) within the uncertainty.
Here, we make the assumption that L bol represents L FIR as young protostars are deeply embedded in gas and dust, and L bol will be dominated by the contribution from L FIR . Hence, if the relation between mass and luminosity is more universal, underestimating or overestimating can respectively underestimate or overestimate the available molecular reservoir.

Comparison with observations
When comparing the derived values with observations, we clearly see that we are underestimating the water emission by at least one to two orders of magnitude (see Fig. 4) and overestimating the SFRs from a factor of a few to two orders of magnitude (see Figs. 3 and 5). We discuss the possible explanations for the difference in SFR in Sect. 4.3 extensively, and here we focus solely on the difference between our estimate and that of Jarugula et al. (2019). The SFR calibration of Jarugula et al. (2019) utilizes the L FIR -SFR relation from Kennicutt & Evans (2012): which, as mentioned in Sect. 4.3, is subject to various uncertainties. This is especially important when considering the IMF in the high-z ULIRGs and HyLIRGs, as found in many studies (e.g., Zhang et al. 2018), adding uncertainty to the calibration. Moreover if we were to apply the calibration from Jarugula et al.
(2019), we would heavily underestimate SFRs towards wellstudied, resolved Galactic clouds, where the relation inferred for water emission and luminosity is ≈3000 times higher than that of Jarugula et al. (2019) (further discussion in Kristensen et al. 2022).
Focusing on the water emission, there are a few factors that could contribute to the observed difference and we discussed some of them in Sect. 4.1. Additionally, one of the reasons for not recovering the emission is that we do not convert 100% of the galactic mass to an emitting source. There is a number of parameters standing in the way, with the star formation efficiency being the most obvious one. Moreover, currently we consider emission only from Class O and Class I protostars. Therefore, when considering emitting components that constitute only a small percentage of a whole galaxy, we are naturally going to lose a certain amount of emission.
In galaxies there are more emitting components than simply protostars. These include photodissociation regions, galactic outflows, and supernovae. Even though their contribution is likely to be lower than that from star formation, their inclusion in calculations is essential in order to fully reproduce the emission, and as such, is a part of planned future improvements.

Conclusions
We extended the galaxy-in-a-box model to relate the predicted molecular emission from forming stars with SFRs. In this paper, we demonstrate the introduced extension and evaluate the derived results for galaxies with L FIR = 10 8 −10 11 L and various levels of star formation activity. We complemented the SFR study by extracting predicted emission for the para-H 2 O 2 02 −1 11 line at 987.927 GHz. Our main results are as follows: -The star formation efficiency and the free-fall-time efficiency have a strong impact on the SFR and emission, whereas the opposite holds for the IMF. -For the most extreme star-forming cases, the galaxy-in-a-box model overestimates the SFRs by up to two orders of magnitude. However, this difference could be lowered depending on the extent to which the current calibrations using L FIR as a star formation tracer underestimate the actual SFR values. -The model underestimates the water emission by up to two orders of magnitude, and especially for the high-z quasars, ULIRGs, and HyLIRGs. -For the moment, the model does not account for additional sources of emission, including supernovae, photodissociation regions, and galactic outflows. Moreover, we need to revisit the derived water emission for Galactic high-massstar-forming regions, as we might miss the bulk of emission. Our estimates deviate from observations and the literature. However, the apparent differences are consistent with expectations in the sense that known sources of emission are not included in the model, and therefore the galaxy-in-a-box model is a promising step toward shedding light on the star-forming properties of galaxies across cosmic time. In the near future, we plan to introduce a number of extensions that will account for other sources and processes that could contribute to the emission. The planned extensions include accounting for galactic outflows -both AGN and starburst driven -, shocks from supernovae, A95, page 9 of 11 and emission from photodissociation regions. Moreover, we are introducing H 2 and high-J CO emission, which is going to be especially important in the JWST era.
To properly account for water emission in our own Galaxy in the future, we will need a new far-infrared probe with the sensitivity of JWST. Such a probe is the planned PRIMA 1 mission. Only then will we be able to fully recover the emission from star-forming clusters in the Galaxy, and properly estimate the contribution from protostars in all stellar mass ranges.