By compiling a comprehensive census of literature studies we investigate the evolution of the main sequence MS of starforming galaxies SFGs in the widest range of redshift 0 lt z lt 6 and stellar mass 10SUP8.5SUP10SUP11.5SUP MSUBSUB ever probed. We convert all observations to a common calibration and find a remarkable consensus on the variation of the MS shape and normalization across cosmic time. The relation exhibits a curvature towards the high stellar masses at all redshifts. The best functional form is governed by two parameters the evolution of the normalization and the turnover mass MSUB0SUBt which both evolve as a power law of the Universe age. The turnover mass determines the MS shape. It marginally evolves with time making the MS slightly steeper towards z  46. At stellar masses below MSUB0SUBt SFGs have a constant specific SFR sSFR while above MSUB0SUBt the sSFR is suppressed. We find that the MS is dominated by central galaxies. This allows to turn MSUB0SUBt into the corresponding host halo mass. This evolves as the halo mass threshold between cold and hot accretion regimes as predicted by the theory of accretion where the central galaxy is fed or starved of cold gas supply respectively. We thus argue that the progressive MS bending as a function of the Universe age is caused by the lower availability of cold gas in haloes entering the hot accretion phase in addition to black hole feedback. We also find qualitatively the same trend in the largest sample of starforming galaxies provided by the IllustrisTNG simulation. Nevertheless we still note large quantitative discrepancies with respect to observations in particular at the highmass end. These can not be easily ascribed to biases or systematics in the observed SFRs and the derived MS.
['galaxies: evolution', 'galaxies: high-redshift', 'galaxies: star formation', 'Astrophysics - Astrophysics of Galaxies']
The main sequence of starforming galaxies across cosmic times
{'P. Popesso, 1 ★ A. Concas, 1 G. Cresci, 2 S. Belli, 3 G. Rodighiero, 4 H. Inami, 5 M. Dickinson, O. Ilbert, 7 M. Pannella, 8 D. Elbaz, 9': "- 1 European Southern Observatory, Karl Schwarzschildstrasse 2, 85748, Garching bei München, Germany\n- 2 INAF-Osservatorio Astronomico di Arcetri, Largo Enrico Fermi 5, 50125, Firenze, Italy\n- 3 Center for Astrophysics Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138\n- 4 Dipartimento di Fisica e Astronomia, Universitá di Padova, Vicolo dell' Osservatorio 3, 35122, Padova, Italy\n- 5 Hiroshima Astrophysical Science Center, 1-3-2 Kagamiyama, Higashi-Hiroshima City, Hiroshima, Japan 739-8511\n- 6 NOIRLab, 950 N. Cherry Ave. Tucson, AZ 85719, USA\n- 7 Laboratoire d'Astrophysique de Marseille, 38 rue Frederic Joliot Curie, 13388 Marseille, France\n- 8 Dipartimento di Fisica e Astronomia, Universitá di Trieste, Via Valerio, 2 - 34127 Trieste, Italy\n- 9 Laboratoire AIM-Paris-Saclay, CEA/DRF/Irfu - CNRS - Université Paris Diderot, CEA-Saclay, 91191 Gif-sur-Yvette, France  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': "By compiling a comprehensive census of literature studies, we investigate the evolution of the Main Sequence (MS) of star-forming galaxies (SFGs) in the widest range of redshift (0 GLYPH<159> 𝑧 GLYPH<159> 6) and stellar mass (10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> ) ever probed. We convert all observations to a common calibration and find a remarkable consensus on the variation of the MS shape and normalization across cosmic time. The relation exhibits a curvature towards the high stellar masses at all redshifts. The best functional form is governed by two parameters: the evolution of the normalization and the turnover mass ( 𝑀 0 ' 𝑡 ' ), which both evolve as a power law of the Universe age. The turn-over mass determines the MS shape. It marginally evolves with time, making the MS slightly steeper towards 𝑧 GLYPH<24> 4 GLYPH<0> 6. At stellar masses below 𝑀 0 ' 𝑡 ' , SFGs have a constant specific SFR (sSFR), while above 𝑀 0 ' 𝑡 ' the sSFR is suppressed. We find that the MS is dominated by central galaxies. This allows to turn 𝑀 0 ' 𝑡 ' into the corresponding host halo mass. This evolves as the halo mass threshold between cold and hot accretion regimes, as predicted by the theory of accretion, where the central galaxy is fed or starved of cold gas supply, respectively. We, thus, argue that the progressive MS bending as a function of the Universe age is caused by the lower availability of cold gas in halos entering the hot accretion phase, in addition to black hole feedback. We also find qualitatively the same trend in the largest sample of star-forming galaxies provided by the IllustrisTNG simulation. Nevertheless, we still note large quantitative discrepancies with respect to observations, in particular at the high mass end. These can not be easily ascribed to biases or systematics in the observed SFRs and the derived MS.  \nKey words: galaxies: evolution - galaxies: star formation - galaxies: high-redshift", '1 INTRODUCTION': 'The Main Sequence (MS) of star forming galaxies (SFGs) is considered one of the most useful tools in modern astrophysics in the field of galaxy evolution. This very tight relation between the galaxy star formation rate (SFR) and the stellar mass (M ★ ) is in place from redshift GLYPH<24> 0 up to GLYPH<24> 6 (Brinchmann et al. 2004; Salim et al. 2007; Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007; Chen et al. 2009; Pannella et al. 2009; Santini et al. 2009; Oliver et al. 2010; Magdis et al. 2010; Rodighiero et al. 2011; Elbaz et al. 2011; Karim et al. 2011; Shim et al. 2011; Whitaker et al. 2012; Zahid et al. 2012; Lee et al. 2012; Reddy et al. 2012; Salmi et al. 2012; Moustakas  \net al. 2013; Kashino et al. 2013; Sobral et al. 2014; Steinhardt et al. 2014; Speagle et al. 2014; Whitaker et al. 2014; Shivaei et al. 2015a; Schreiber et al. 2015; Tasca et al. 2015; Lee et al. 2015; Kurczynski et al. 2016; Erfanianfar et al. 2016; Santini et al. 2017; Pearson et al. 2018; Belfiore et al. 2018; Popesso et al. 2019a,b; Leslie et al. 2020; Thorne et al. 2021; Leja et al. 2022; Daddi et al. 2022). The evolution of its normalization, slope and scatter have been largely studied in the past decade. It is now well established that the normalization declines significantly but smoothly as a function of redshift, likely on mass-dependent timescales (see also Speagle et al. 2014), rather than being driven by stochastic events like major mergers and starbursts (Oemler et al. 2017). More uncertain is the precise redshift dependence of such evolution, which is often expressed as / \' 1 , 𝑧 \' 𝛾 , with 𝛾 varying from 1.9 to 3.7 (Speagle et al. 2014;  \nWhitaker et al. 2014; Schreiber et al. 2015; Ilbert et al. 2013; Pearson et al. 2018; Leslie et al. 2020; Thorne et al. 2021; Leja et al. 2022). This is mainly due to the uncertainty in deriving the evolution of the exact shape of the relation, which is still matter of intense debate. Several studies point to a power law shape, SFR / M 𝛼 ★ , both in the local Universe (Peng et al. 2010; Renzini & Peng 2015) and at high redshift (Speagle et al. 2014; Rodighiero et al. 2014; Kurczynski et al. 2016; Pearson et al. 2018). Other works suggest that the relation exhibits a curvature towards the high mass at low (Popesso et al. 2019b) and high redshift (Whitaker et al. 2014; Schreiber et al. 2015; Lee et al. 2015; Tasca et al. 2015; Tomczak et al. 2016; Popesso et al. 2019a; Leslie et al. 2020; Thorne et al. 2021; Leja et al. 2022). Also the scatter around the relation is quite debated with very conflicting results from the literature. Some works report a quite constant scatter of 0.2-0.3 dex from low to moderately high masses (e.g. Schreiber et al. 2015; Noeske et al. 2007; Elbaz et al. 2007). Others, instead, report a decrease of the scatter from very low (10 8 𝑀 GLYPH<12> ) to moderate (10 10 𝑀 GLYPH<12> ) stellar masses at different redshifts (Willett et al. 2015; Santini et al. 2017; Boogaard et al. 2018). Few others observe an increase as a function of the stellar mass from 0.3 dex at 10 10 to 0.5 dex at 10 11 GLYPH<149> 5 𝑀 GLYPH<12> , from low to high redshift(e.g. Guo et al. 2013; Popesso et al. 2019b,a; Sherman et al. 2021).  \nMost of this discrepancy is ascribable to how SFGs are selected in the first place, on the SFR estimators and on the method of localization of the MS. Popesso et al. (2019a) show that color-color selection of SFGs leads to the exclusion of red dusty star-forming galaxies in particular at 𝑧 ¡ 1 GLYPH<149> 5 and to a steep MS. Tomczak et al. (2016) show that without any selection, the MS is bending more significantly in the local than in the distant Universe (see also Katsianis et al. 2016). In addition, the use of different SFR indicators might lead to systematic biases. Part of the UV emission originating from the young star population is absorbed by dust and re-processed at infrared wavelengths. Such emission alone can provide a measure of the SFR only if corrected for this absorption. However, the measure of the dust attenuation is still uncertain because of the degeneracy between age and reddening, the assumption on galaxy metallicity and SF histories, and the parametrization of the extinction curve (e.g. Meurer et al. 1999; Dale et al. 2009; Dunlop et al. 2017; Bourne et al. 2017). With the launch of the Spitzer and Herschel satellites in 2003 and 2009, respectively, it became possible to measure the mid and far-infrared emission for statistical samples of galaxies up to 𝑧 GLYPH<24> 3 in the most studied deep fields (e.g. Lutz 2014). Such measurements allowed finally combining the unobscured UV emission and the reprocessed far-infrared (FIR) component and calibrating them as a star formation rate indicator. Nevertheless, also such an approach might lead to systematics and uncertainty. For instance, the contamination by AGN and the overestimation of the SFRfor starburst galaxies make of the mid-infrared emission, based e.g. on Spitzer MIPS 24 𝜇 m data, a less accurate indicator with respect to longer wavelengths (Elbaz et al. 2011; Nordon et al. 2010). Nonetheless, such contamination can affect, to some extent, also the emission in the far-infrared, if the nuclear activity was particularly high in the last 100 𝑀𝑦𝑟 of a galaxy\'s life. More recently, the accuracy of the SFR based on the combination of UV and IR components has been questioned by the results based on sophisticated SED fitting codes that try to reconstruct the galaxy star formation histories (e.g. Thorne et al. 2021; Leja et al. 2022). Navigating through all these differences makes it very difficult to understand if the literature has reached an overall consensus on the shape and evolution of the Main Sequence across cosmic times. Furthermore,  \nthe magnitude of these effects precludes robust interpretations of derived MS properties.  \nTo overcome this problem, and so constrain the MS evolution and systematic errors, Speagle et al. (2014, hereafter S14) have compiled 64 MS observations from 25 studies published in the period 2007-2014, spanning 𝑧 GLYPH<24> 0 - 6, and converted them to the same absolute calibration. These MS estimates have been taken from a variety of fields, selected using different methodologies, including both stacked and non-stacked data, and estimated with a variety of SFR indicators. By calibrating consistently all datasets, S14 determine the MS best fit as a power law, with a slope marginally evolving with redshift. The purpose of such an experiment is not to provide the "true" Main sequence. Indeed, there might still be biases and limitations related to the observational techniques, that this approach is not capable of erasing or correcting (see for instance the discussion in Katsianis et al. 2020). Rather, it aims at understanding whether the MS estimates, reported in the literature, lead to a substantial consensus once they are brought to the same calibration and all systematics are considered. It offers also the advantage that the resulting MS can be easily re-calibrated, if necessary.  \nIn this paper, we adopt the same approach of S14, to bring all estimates of the literature to a common calibration to check whether the community has reached or is far from reaching such a consensus. To this aim, we extend the collection of MS determination of S14 to the most recent estimates of the MS based on a variety of SFR indicators (from the combination of Spitzer mid-infrared and Herschel far-infrared data and UV emission, to SED fitting techniques) and methodology. These include additional 27 publications from 2014 to 2022, for a total of GLYPH<24> 120 MS measurements at 0 GLYPH<159> 𝑧 GLYPH<159> 6, that we convert to the same IMF and calibrate in a consistent framework. We collect a sample of about GLYPH<24> 1500 consistently calibrated data points, that express the location of the MS as a function of stellar mass and time, and check whether there is consensus in the evolution of the MS across cosmic time.  \nThe paper is structured as follows. Section 1 describes the MScompilation. Section 2 presents our best-fitting procedure. Section 3 shows our results. Section 4 provides a comparison of our results with previus findings and throretical predictions, while Section 5 lists our conclusions. We assume a Λ CDM cosmology with Ω 𝑀 = 0 GLYPH<149> 3, Ω Λ = 0 GLYPH<149> 7 and 𝐻 0 = 70 km/s/Mpc , and a Kroupa IMF throughout the paper.', '2 THE MS COLLECTION': 'In this work, we focus on MS estimates that have been published after 2014. To consider in the analysis the MS measurements that have been published in the period 2007-2014, we include here the MSrelation determined by S14 as a result of a compilation of 64 MS estimates collected in 25 papers (see Table 3 of S14). We consider, in particular, the fit n. 64, which S14 provide as their reference MS. To this, we add all the MS relations that satisfy the following criteria:  \n- (i) Includes a published M ★ - SFR or M ★ - sSFR (sSFR=SFR/M ★ ) relation (slope 𝛼 and normalization 𝛽 ) or otherwise analogous quantities;\n- (ii) Fit(s) includes more than two data points (if stacked) or 50 galaxies (if directly observed). This is required to avoid biases resulting from small number statistics;\n- (iii) Stacked points must provide mean or median of more than 25 points, to avoid large uncertainties due to low number statistics; (iv) Published after 2014.  \nTable 1. Col. 1: Reference. Col. 2: Assumed stellar initial mass function. Col. 3: Star formation rate indicator. Col. 4: redshift range analyzed Col. 5: Assumed cosmology. Col. 6: Selection methods used for the parent samples (see Appendix or individual papers for more details). Col. 7: type of data retrieved in the paper: stacked stands for stacked SFR data-points as a function of stellar mass and redshift, data-points for average or median SFR data-points as a function of stellar mass and redshift, and best fit stands for best-fit parameters expressing the MS shape as a function of stellar mass over the retrieved stellar mass and redshift ranges. Col. 8: Extinction curve reference: the details of the extinction curve implemented in MAGPHYS can be found in da Cunha et al. (2008), C00 refers to Calzetti et al. (2000), CF00 to Charlot & Fall (2000), KC13 to Kriek & Conroy (2013) and N09 to Noll et al. (2009). Col 9: indicates whether the data have been included in ( 3 symbol) or excluded from ( GLYPH<0> symbol) the current analysis. The reasons for excluding an individual dataset from the analysis are given in Section 2. and In Appendix ?? with more details.  \nThe 27 publications that are retrieved considering these criteria, in addition to S14, are listed in Table 1, together with the used IMF, SFR indicator, redshift range, cosmological parameters, SFGs selection method, and extinction curve. In Appendix, we provide a full description of the individual publication data and of their calibration.  \nAll the considered MS estimates are based on the deepest UV, optical, IR and radio surveys ever realized on the CANDELS, COSMOS and ECDFS fields at intermediate and high redshift, and on the WISE, GALEX and optical spectroscopy in the SDSS area at 𝑧 GLYPH<24> 0. Of these, 10 publications (Heinis et al. 2014; Whitaker et al. 2014; Chang et al. 2015; Lee et al. 2015; Ilbert et al. 2015; Schreiber et al. 2015; Tomczak et al. 2016; Pearson et al. 2018; Lee et al. 2018; Popesso et al. 2019a) out of 27 are based on an SFR indicator given by the combination of UV and IR data. Additional 3 are based only on far-infrared PACS data (Rodighiero et al. 2014; Erfanianfar et al. 2016; Popesso et al. 2019a). Other 3 (Renzini & Peng 2015; Belfiore et al. 2018) are based on 𝐻𝛼 derived SFR taken from the SDSS spectroscopic dataset (Brinchmann et al. 2004) or from 3D-HST data at higher redshift (de los Reyes et al. 2015). Additional 9 include SFR derived through SED fitting technique based on different fitting methods and star formation history recon-  \ns (Tasca et al. 2015; Salmon et al. 2015; Santini et al. 2017; Kurczynski et al. 2016; Iyer et al. 2018; Barro et al. 2019; Thorne et al. 2021; Sherman et al. 2021; Leja et al. 2022). Of the rest, one is based on the evolution of the galaxy stellar mass function of starforming galaxies (Davidzon et al. 2018), and one is based on deep radio data (Leslie et al. 2020). S14 is based on a collection of different SFR estimators. According to the definition of S14, all of the considered publications but two are based on "mixed" methods (see next section) for the selection of SFGs. These include color-color techniques (BzK, UVJ, and NUVRJ), 2 𝜎 clipping, and bimodality in the SFR-M ★ plane. Only 3 publications have an MS obtained for blue and active galaxies selected in the UV or with the Lyman Break Galaxy (LBG) technique.  \nDifferently from S14, we include in our study also the MS estimate at 𝑧 GLYPH<24> 0. S14 discuss the inability to distinguish a "best" MS fit among the available 𝑧 GLYPH<24> 0 estimates obtained before 2014. Thus, they decide not to include the local MS estimates. Popesso et al. (2019b) discuss extensively all the datasets available in the local Universe, selection effects and different SFR estimators and the level of agreement between the different estimates. On the basis of these results, we conclude that the inter-publication scatter of the local MS estimates available in the period 2014-2022 is comparable to that  \nfound at higher redshifts. Indeed, no different level of agreement is found as a function of redshift, once all MS estimates are brought to the same calibration (see next section for more details).  \nFrom each publication either we take the mean or median SFR data points or staked SFR data points at the observed stellar mass and redshift without any extrapolation or interpolation, or, if these are not available, we use the provided MS best-fit parameters at a given redshift to estimate the MS in the provided stellar mass range in bins of 0.15 dex in stellar mass. The bin width is chosen to be representative of the average stellar mass error of the considered papers. The mass ranges have either been taken directly from the paper in question or estimated based on the data included in the relevant fits, rounded to the nearest 0 GLYPH<149> 1 dex (see Appendix for a detailed description of the data taken from each publication). This leads to GLYPH<24> 120 determinations of the MS for a sample of GLYPH<24> 1500 data points of MS SFR as a function of stellar mass and redshift or time ( 𝑆𝐹𝑅 \' 𝑀 GLYPH<3> GLYPH<150> 𝑧 \' or 𝑆𝐹𝑅 \' 𝑀 GLYPH<3> GLYPH<150> 𝑡 \' ). The data collected here encompass the widest range in redshift (0 GLYPH<159> 𝑧 GLYPH<159> 6), stellar mass (10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> ), and SFR (0 GLYPH<149> 01 GLYPH<0> 500 𝑀 GLYPH<12> yr GLYPH<0> 1 ) available in the literature, and intend to give a census of most of the techniques and methods used to derived the MS location. Following the example of S14, we present what we hope is the broadest and most accurate census of MS observations to date.', '2.1 MS calibration': 'As underlined in S14, several aspects need to be taken care of when comparing different MS estimates and before attempting any analysis:  \n- (i) Initial Mass Function;\n- (ii) SFR estimator;\n- (iii) SPS model;\n- (iv) cosmology;\n- (v) emission line effect in the estimate of SFR and M ★ ;\n- (vi) Star formation histories (SFH);\n- (vii) dust extinction curve;\n- (viii) photo-z biases;\n- (ix) SED fitting procedures;\n- (x) selection effects due to different SFG population selection.  \nA detailed discussion in S14 points out that only the points i) , ii) and iii) lead to relevant corrections when calibrating all the MS estimates to a common ground. Instead, different cosmologies ( iv ) have relatively negligible effects ( GLYPH<159> 0 GLYPH<149> 05 dex) at 0 GLYPH<159> 𝑧 GLYPH<159> 6, which is the same redshift range considered here. The effect of emission lines ( v ) in the estimate of M ★ and SFR is relevant only when the multiwavelength information are limited to few photometric bands. They are, instead, negligible for datasets like CANDELS, and COSMOS, as those considered here.  \nSimilarly to S14, we choose not to adjust our results for differences in assumed star formation history ( vi ), dust attenuation curves ( vii ), possible photo-z biases ( viii ), or differences in SED fitting procedures ( ix ). This might differ substantially (see the extensive discussion in Simha et al. 2014; Acquaviva et al. 2015; Salmon et al. 2015; Ciesla et al. 2017; Iyer & Gawiser 2017; Theios et al. 2019; Carnall et al. 2019; Lower et al. 2020; Thorne et al. 2021; CurtisLake et al. 2021). Testing and correcting the effect due to such differences would require checking the galaxy SFR catalogs used to build the MS of the individual publications. In none of the considered cases, such information is publicly available. Only the stacked, averaged, or fitted MS is provided, which prevents performing any correction. In principle, the effect of the different assumptions and  \nimplementations could be predictable. There is a large recent literature of such attempts based on the combination of dust radiative transfer models and simulations (Baes et al. 2020; Trayford et al. 2020; Lower et al. 2020; Narayanan et al. 2021; Lovell et al. 2021; Katsianis et al. 2021). However, such corrections might introduce biases dependent on the different assumptions of the models. Thus, we choose not to apply such corrections. As shown at the end of the next section, neglecting the possible effects of the mentioned differences does not lead to an increase of the resulting inter-publication scatter.  \nWe apply a correction only for 2 publications (Thorne et al. 2021; Leja et al. 2022), for which the authors explicitly indicate a systematic difference in the estimate of the stellar mass with respect to previous works of our collection. This correction is applied not because we consider these quantities wrong or inaccurate, but only to bring all MS estimates to the same stellar mass and SFR scale (see next paragraph and Appendix A). All other publications based on the SED fitting technique, do not report systematic differences in their stellar mass and SFR estimates. Thus, no further corrections are applied.  \nSFG selection methods ( x ) can lead to substantially different MS slopes. S14 distinguish between "bluer", "mixed", and "nonselective" techniques: "bluer" methods are based on a simple color cut to select blue galaxies, "mixed" ones involve a more sophisticated distinction between star-forming and quiescent galaxies, while "non-selective" ones do not apply any selection. In particular, S14 find that "bluer"-based ("non-selective"-based) MS slopes are biased towards values closer to unity (zero), with respect to "mixed"based slopes.  \nTo bring every relation onto a common framework, we apply the following steps. We use the Equations reported in Section 2.2 to convert the M ★ values to a Kroupa IMF, and the SFR estimates to the Kennicutt & Evans (2012, hereafter KE12) calibration (based on the Kroupa IMF). The choice of the KE12 calibration is dictated by the fact that it is the only one available in the literature that calibrates consistently all the SFR indicators considered here, with the exclusion of the SED fitting technique. Thus, it is the obvious choice to bring most of the SFR estimates to a common framework. The available SFRs based on the SED fitting technique included in the collection, are compared case by case with the KE12 calibration to check for consistency, as described in Appendix A. The KE12 calibration is based on data of star forming regions in a nearby galaxy (Murphy et al. 2011; Hao et al. 2011). To check that the SFR indicators calibrated according to KE12 provide consistent estimates also at higher redshift, we perform several tests in Appendix B. The result of such tests is that SFRs derived through 𝐼 𝑅 , 𝑈𝑉 , 𝐼 𝑅 , 𝐻𝛼 and radio emission at 1.4 GHz according to the KE12 calibration (see next section), are all consistent up to 𝑧 GLYPH<24> 3 with a scatter that varies from 0.2 to 0.25 dex.  \nWe adjust for differences in cosmology using the \' ℎGLYPH<150> Ω 𝑀 GLYPH<150> Ω Λ \' = \' 0 GLYPH<149> 7 GLYPH<150> 0 GLYPH<149> 3 GLYPH<150> 0 GLYPH<149> 7 \' WMAP concordance cosmology (Spergel et al. 2003). Furthermore, no additional dust correction is required, as the dust has been corrected for in the considered publications, and different extinction curves do not have a significant impact on the MS determination (see S14). Finally, as only 3 (Tasca et al. 2015; Heinis et al. 2014; Thorne et al. 2021) out of the 27 MS estimates considered here are consistent with a "bluer" selection, we do not correct in this respect, but we discuss possible biases on the result.  \nAs pointed out in Popesso et al. (2019b), another source of small discrepancy among different MS estimates is the method used to determine the MS location, e.g. as the mean or median  \nSFR of the SFG population in the MS region. A correction can be made in this respect under the assumption that the SFR distribution at fixed M ★ is log-normal in the MS region. This assumption is justified by several works in literature, which find a log-normal distribution at any stellar mass (Daddi et al. 2007; Rodighiero et al. 2011; Schreiber et al. 2015; Popesso et al. 2019b,a; Leja et al. 2022). In this case, the peak of the distribution coincides with the median SFR. The mean SFR is always larger than the median by an offset that depends only on the dispersion of the distribution. For a dispersion of GLYPH<24> 0 GLYPH<149> 3 dex, the correction is GLYPH<24> 0 GLYPH<149> 1 dex and it increases to GLYPH<24> 0 GLYPH<149> 17 dex for a dispersion of GLYPH<24> 0 GLYPH<149> 5 dex. As pointed out in the Introduction, many discrepant results in the literature do not allow us to reach a consensus on the value of the MS scatter as a function of the stellar mass. Thus, we assume an average dispersion of 0.3 dex and apply the corresponding correction in the further analysis to all MS estimates based on the median SFR to convert them into mean values.', '2.2 The main corrections': 'We describe here the corrections applied to bring the different MS estimates to a common framework:', 'GLYPH<15> IMF correction': 'We apply the IMF offsets to stellar masses as in S14 with the form:  \n𝑀 ★GLYPH<150> 𝐾 = 1 GLYPH<149> 06 𝑀 ★GLYPH<150> 𝐶 = 0 GLYPH<149> 62 𝑀 ★GLYPH<150> 𝑆 GLYPH<150> (1)  \nwith the subscripts referring to Kroupa, Chabrier, and Salpeter IMFs, respectively. These correspond to stellar mass offsets of 0 GLYPH<149> 03 and GLYPH<0> 0 GLYPH<149> 21 dex, for the Chabrier and the Salpeter IMF, respectively.', 'GLYPH<15> SFR indicator correction': "All SFR estimates have been converted to the calibration of KE12, which is based on the Kroupa IMF. The KE12 prescriptions are different depending on whether 𝜈𝐿 𝜈 is estimated in the FUV ( GLYPH<24> 1300 GLYPH<0> 1700Å) or NUV ( GLYPH<24> 2300 GLYPH<0> 2800Å). Following Tables 1 and 2 of KE12, we make use of the following calibrations:  \n𝑆𝐹𝑅 ' 𝐹𝑈𝑉 , 𝐼 𝑅 ' = 1 GLYPH<149> 71 GLYPH<1> 10 GLYPH<0> 10 ' 𝐿 𝐹𝑈𝑉 , 0 GLYPH<149> 46 GLYPH<3> 𝐿 𝐼 𝑅 ' (2)  \n𝑆𝐹𝑅 ' 𝑁𝑈𝑉 , 𝐼 𝑅 ' = 2 GLYPH<149> 06 GLYPH<1> 10 GLYPH<0> 10 ' 𝐿 𝑁𝑈𝑉 , 0 GLYPH<149> 27 GLYPH<3> 𝐿 𝐼 𝑅 ' (3)  \n𝑆𝐹𝑅 ' 𝐻𝛼 ' = 2 GLYPH<149> 06 GLYPH<1> 10 GLYPH<0> 8 ' 𝐿 𝐻𝛼 ' (4)  \n𝑆𝐹𝑅 ' 𝐼 𝑅 ' = 1 GLYPH<149> 49 GLYPH<1> 10 GLYPH<0> 10 ' 𝐿 𝐼 𝑅 ' (5)  \n𝑆𝐹𝑅 ' 1 GLYPH<149> 4 𝐺𝐻𝑧 ' = 2 GLYPH<149> 43 GLYPH<1> 10 5 ' 𝐿 1 GLYPH<149> 4 𝐺𝐻𝑧 ' (6)  \nwhere 𝐿 𝐼 𝑅 , 𝐿 𝐹𝑈𝑉 , 𝐿 𝑁𝑈𝑉 and 𝐿 𝐻𝛼 are the luminosities (in solar luminosity units) estimated in the IR (range 8 GLYPH<0> 1000 𝜇 m), FUV and NUV, 𝐻𝛼 emission line and radio emission (1.4 GHz), respectively. The SFRs based on SED modeling, instead, are corrected only for the IMF, according to the derivative of Eq. 1. The correction to a common SFR indicator varies from 0.05 to GLYPH<0> 0 GLYPH<149> 2 dex. After applying the calibration, the local MS of Chang et al. (2015) is still systematically lower with respect to the other local MS estimates by 0.15 dex (see Popesso et al. 2019b for a detailed comparison with UV , IR based and H 𝛼 based SFR). Thus, we correct for this  \noffset as indicated in Popesso et al. (2019b) before including the MS estimate of Chang et al. (2015) into our sample.  \nAs already pointed out in Popesso et al. (2019b), Pearson et al. (2018) report a systematic offset of their SFR of GLYPH<24> 0.4 dex below all other SFR estimates considered here at the same redshift. As shown in Popesso et al. (2019b), their MS lies below all other determinations at more than 1 𝜎 at all redshift. Elbaz et al. (2011) show that SPIRE and PACS SFR estimators lead to consistent results. So we conclude that the discrepancy must be related to the deblending technique of the SPIRE detections and the SED fitting technique applied in Pearson et al. (2018, see their Appendix C for an extensive discussion). Since the problem might be related to an over-deblending issue rather than to the SFR indicator, we decide not to correct the SFR of Pearson et al. (2018) and to exclude those MS estimates from the dataset.  \nIn addition, two of the considered publications (Thorne et al. 2021; Leja et al. 2022) based on two different SED fitting algorithms, ProSpect (Robotham et al. 2020), and Prospector (Leja et al. 2017), report lower values of SFR and larger values of stellar masses with respect to the rest of the MS estimates considered here. This discrepancy is due to a different reconstruction of the galaxy star formation history with respect to other SED fitting methods. The main difference consists in the fact that both codes include the contribution of an extra component of stars older than 100 𝑀𝑦𝑟 , which would affect both the stellar mass and the SFR estimates. The effect on the stellar masses is obvious and it is reported to be 0.2 dex for Thorne et al. (2021) and 0.3 dex for Leja et al. (2022). Leja et al. (2022), in particular, report also that the derived SFRs tend to be lower at fixed stellar mass with respect to 𝐼 𝑅 or 𝑈𝑉 , 𝐼 𝑅 derived SFR in the redshift range 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3. They ascribe such a difference to the contribution of the extra component of old stars to the dust heating, which would increase the galaxy IR emission. They also find that such contribution is redshift dependent, being larger at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3 than at lower and higher redshifts. Leja et al. (2022) conclude that 𝐼 𝑅 or 𝑈𝑉 , 𝐼 𝑅 derived SFRs are overestimated at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3, because not all the IR emission is due to dust heating by young stars and, thus, to star formation.  \nThe KE12 calibration does take into account that only a fraction of the 𝐼 𝑅 emission is due to star formation. However, it assumes that such a fraction is the same at all redshifts. If this assumption is not valid, it implies that at higher redshift, the 𝐼 𝑅 or 𝐼 𝑅 , 𝑈𝑉 SFR indicators based on the KE12 calibration might lead to over or underestimated SFRs with respect to the radio or H 𝛼 derived SFRs, based on the same calibration. This is because the contribution of an extra old star component might significantly affect the IR emission, but it has rather a negligible effect on the H 𝛼 and radio emission, which are dominated by HII region nebular emission and non-thermal (synchrotron) emission due to Type II and Type Ia Supernovae, respectively. As shown in the Appendix B, we do not observe such discrepancy at any redshift. Instead, the KE12 calibrated SFR derived through the radio, H 𝛼 and 𝐼 𝑅 , 𝑈𝑉 emission are consistent from 𝑧 GLYPH<24> 0 up to 𝑧 GLYPH<24> 3. Thus, we do not find clear evidence of an evolution of the fraction of 𝐼 𝑅 emission contributing to the SFR. In addition, we point out that Thorne et al. (2021) show that once the stellar mass discrepancy is taken into account, the SFRs derived with ProSpect are perfectly consistent with the radio emission based SFRs of Leslie et al. (2020). We perform the same exercise between the Prospector MS estimates of Leja et al. (2022) and those of Whitaker et al. (2014) based on 𝐼 𝑅 , 𝑈𝑉 emission on the same 3DHST galaxy sample in the CANDELS fields (see Fig. A1 in Appendix A). Once the stellar masses are corrected for the discrepancy (0.3 dex), the MS estimates are consistent within  \n<!-- image -->  \nFigure 1. Left panel: logSFR versus Universe age in several bins of stellar masses. The data points indicate the SFR based on the MS estimates collected in this work. The solid lines show the best linear fit as in Eq. 9. Data points and lines are color-coded as a function of the stellar mass bin as indicated in the figure. For clarity, the relations are artificially displaced by 0.4 dex from one another. Right panel: 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' (upper panel) and 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' (lower panel) as a function of 𝑙𝑜𝑔𝑀 ★ . The red solid lines in both panels indicate the best fit relations of Eq. 9. The dashed line in the bottom panel shows the liner fit approximation as proposed in S14.  \n<!-- image -->  \n- 1-1.5 𝜎 at all stellar masses. Thus, we include the Thorne et al. (2021) and Leja et al. (2022) MS estimates by correcting the stellar masses by the reported discrepancy, and we leave unaltered the SFR estimates.\n- GLYPH<15> Cosmology correction : this is calculated as the ratios between luminosity distance, 𝑑 𝐿 ' 𝑧 ' , derived from two different cosmologies, and, given the observed redshift range of a sample, applying a 𝑑 2 𝐿 correction at the expected median 𝑧 of galaxies in the sample. S14 estimate also first-order volume effects. However, they find that this account for a negligible effect in all cases. Thus, we do not take it into account.  \nAfter applying the calibration, we obtain an inter-publication scatter of 0.08 dex per bin of time and stellar mass. We observe also that the scatter is not much dependent on the SFR indicator or methodology used in the different publications. If we limit the analysis only to publications based on 𝐼 𝑅 or 𝑈𝑉 , 𝐼 𝑅 based SFRs or only to those based on the SED fitting technique, we do not observe any scatter variation.", '2.3 Selection effects': 'While efficient at selecting SFGs, most selection techniques differ from each other and do not all select the same population. S14 discuss extensively that the selection of blue objects preferentially select actively star-forming, non-dusty galaxies and exclude a large percentage of galaxies that are classified as SFGs via other selection mechanisms (e.g. color-color selection). This leads on average to larger MS slopes than the ones retrieved with other selection methods. One of the methods to retrieve blue sources is the Lyman break technique (Steidel et al. 1999; Stark et al. 2009; Bouwens et al. 2011), used to select high𝑧 Lyman-break galaxies (LBGs). In the list of publications considered here only Tasca et al. (2015)  \napply the LBGs selection technique. Tomczak et al. (2016) compare the results of Tasca et al. (2015) with those based on the "mixed" selection method (see also discussion below) and find very good agreement. We conclude that the results of Tasca et al. (2015) seem less affected by the bias of the bluer selection discussed by S14. Similarly, Heinis et al. (2014) apply a UV selection to identify distant SFGs. Nevertheless, their estimates are not scattering significantly with respect to the other relations. Thus, we conclude that also in this case, the UV selection does not affect significantly the slope of the relation. Thorne et al. (2021) applies a sSFR cut, which might similarly select blue galaxies. However, the authors report a very good agreement with the MS estimates of Leslie et al. (2020), which is based on the radio selection. Thus, we conclude that also in this case the bias is negligible.  \nAll other MS estimates included in our analysis are based on methods that S14 classify as "mixed". Such techniques include redder objects in the selection, thus considering also a large portion of the SFG population dominated by dust. These methods are considered to provide a more physical distinction between SFGs and quiescent galaxies (Ilbert et al. 2013; Schreiber et al. 2015). Among these "mixed" methods, we include the color-color selection based on the rest-frame (U-V) GLYPH<0> (V-J) and (M 𝑁𝑈𝑉 -M 𝑅 ) GLYPH<0> (M 𝑅 -M 𝐽 ) absolute colors, the 2 𝜎 clipping, the bimodality between SFGs and quiescent galaxies in the SFR-M ★ plane. While these different selection methods do not seem to affect the average observed SFRs across different publications, as pointed out by S14, they do seem to influence the derived slopes and the intrinsic scatter of the MS. Popesso et al. (2019a) show that all these methods agree very well over most of the stellar mass range considered. However, small discrepancies can be observed due to little selection biases at very high stellar masses towards high redshift. Namely, selections as the (U-V) GLYPH<0> (V-J) color selection tend to exclude part of the high mass star-forming galaxies at relatively lower SFR with respect to the MS  \nFigure 2. The upper panel shows the MS shape as a function of the Universe age from 1 to 13 Gyrs. The data points show the MS estimates collected and calibrated in this work. The solid lines indicate the best fit relation obtained with Eq. 10. Points and lines are color-coded as a function of time, as indicated in the panel. For clarity, the MS relations are artificially displaced by 0.4 dex from one another. The bottom panel shows the residual distribution as a function of stellar mass. The dashed red line shows the 0 level, corresponding to the best fit value. The red points indicate the mean value of the residuals, while the dotted and dashed-dotted lines indicate the 1 and 3 𝜎 regions, respectively.  \n<!-- image -->  \nlocation. This would lead to a steeper MS at higher redshift with respect to the low redshift relation and thus a more significant evolution. The effect is an offset of GLYPH<24> 0.15 GLYPH<0> 0.2 dex at M ★ of 10 11 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> . Instead, the selection based on the (M 𝑁𝑈𝑉 -M 𝑅 ) GLYPH<0> (M 𝑅 -M 𝐽 ) absolute colors is less prone to this selection bias, as it allows to select all galaxies populating the log-normal distribution around the MS location (Popesso et al. 2019a; Ilbert et al. 2015).  \nWe point out that the MS estimates affected by this bias are the stacked points based on the CANDELS field dataset (Whitaker et al. 2014; Schreiber et al. 2015; Tomczak et al. 2016). Due to the very small volume sampled by the CANDELS fields, these include less than 15 galaxies per stacked point at stellar masses larger than 10 11 𝑀 GLYPH<12> . Since this is below the limit required in Section 2, those stacked points are anyhow not included in our analysis.', '3 FITTING THE MAIN SEQUENCE': 'In this work, we adopt two approaches to fit the MS relation. The first one follows S14 and consists in looking for the functional form that best describes the MS relation. As an alternative approach, we investigate the functional form proposed by Lee et al. (2015), which has the advantage of being expressed with parameters that have a physical description. In the following procedure, after correcting stellar masses and SFRs with the calibration described above, we take into account the redshift and mass ranges of each study, and we include in the fitting procedure only objects or stacked points actually observed at given mass and redshift without any extrapolation  \nor interpolation. These mass ranges have either been taken directly from the paper in question or estimated based on the data included in the relevant fits, rounded to the nearest 0 GLYPH<149> 1 dex after excluding outlying points (see Appendix A for a detailed description of the data taken from each publication).', '3.1 S14 approach': "In this Section, we describe how we fit the MS and retrieve the functional form that best expresses the evolution of the slope and normalization as a function of time. We follow a revised version of the S14 approach. We proceed to fit the evolution of the SFR provided by each individual MS as a function of time:  \n𝑙𝑜𝑔𝑆𝐹𝑅 ' 𝑡 ' = 𝛼 𝑖 𝑡 , 𝛽 𝑖 GLYPH<150> (7)  \nThe choice of fitting as a function of time rather than redshift is mainly practical, as a straightforward linear fit works very well for the time variable, while a more complicated functional form is required as a function of redshift.  \nBy fitting 𝛼 𝑖 's and 𝛽 𝑖 's for a grid of M ★ , as shown in Fig. 1, we can derive a function of the form:  \n𝑙𝑜𝑔𝑆𝐹𝑅 ' 𝑡GLYPH<150> 𝑙𝑜𝑔𝑀 ★ ' = 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' 𝑡 , 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' GLYPH<150> (8)  \nassuming a given parametrization for 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' and 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' .  \nAs shown in the right panel of Fig. 1, the slope 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' depends linearly on 𝑙𝑜𝑔𝑀 ★ , consistently with S14. Conversely, the best fitting form for 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' is not a simple linear dependence but it requires a quadratic form. This is because, differently for nearly all the MS compiled by S14, most of the MS estimates included here find a bending of the MS towards large stellar masses. Thus, the best fitting functions for 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' and 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' are, respectively:  \n𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' = 𝑎 , 𝑎 𝑙𝑜𝑔𝑀  \n𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' = 𝑏 , 𝑏 𝑙𝑜𝑔𝑀 ★ , 𝑏 𝑙𝑜𝑔 𝑀 ★  \n0 1 ★ 0 1 2 2 (9)  \nwhich gives:  \n𝑙𝑜𝑔𝑆𝐹𝑅 ' 𝑡GLYPH<150> 𝑙𝑜𝑔𝑀 ★ ' = ' 𝑎 𝑡 , 𝑏 ' 𝑙𝑜𝑔𝑀  \n, 𝑏 2 𝑙𝑜𝑔 𝑀 ★ , ' 𝑏 0 , 𝑎 0 𝑡 ' GLYPH<149> (10)  \n1 1 ★ 2  \nEq. 10 differs from the best fitting function of S14 for the quadratic term 𝑏 2 𝑙𝑜𝑔 2 𝑀 ★ , which is time-independent. The slope of the MS is driven by the combination of the quadratic term and the linear term that sets the faint end slope, ' 𝑎 1 𝑡 , 𝑏 1 ' 𝑙𝑜𝑔𝑀 ★ , which evolves with time. The normalization of the relation depends linearly on time through the term ' 𝑏 0 , 𝑎 0 𝑡 ' .  \nWe use the functional form retrieved in Eq. 10 to fit the whole MS dataset without binning in stellar masses and time. We limit, however, the fitting procedure to the 10 8 GLYPH<149> 7 GLYPH<0> 10 11 GLYPH<149> 3 𝑀 GLYPH<12> range, where the MS included in our collection were declared to have high completeness in stellar mass. The best fit parameters of Eq. 9 for 𝛼 ' 𝑙𝑜𝑔𝑀 ★ ' and 𝛽 ' 𝑙𝑜𝑔𝑀 ★ ' are used as first guess for the fitting procedure. The final best-fit parameters are given in Table 2 and the best-fit MS is shown as a function of time in Fig. 2. The bottom panel shows the residual distribution of the MS estimates with respect to the best fit at a given time and stellar mass. The scatter around the best fit is 0.09 dex. The MS is found to bend towards large stellar masses at all times due to the time-independent quadratic term. The time variation of the faint end slope makes the relation steeper at early epochs. Nevertheless, the evolution of the slope is marginal after the first GLYPH<24> 4 Gyrs, consistently with the results of S14, and as found in Popesso et al. (2019a). The normalization of the relation", '8': "<!-- image -->  \nFigure 3. Left panel: evolution of 𝛾 (upper panel), 𝑙𝑜𝑔 ' 𝑀 0) (central panel) and 𝑙𝑜𝑔 ' 𝑆𝐹𝑅 𝑚𝑎𝑥 ' (bottom panel) as a function of time. The red solid lines in each panel indicate the best fit relations: a constant equal to 1 for 𝛾 and the best fit of Eq. 12 and Eq. 13 for 𝑀 0 and 𝑆𝐹𝑅 𝑚𝑎𝑥 , respectively. Right panel: the upper panel shows the MS shape as a function of the Universe age from 1 to 13 Gyrs. The data points show the MS estimates collected and calibrated in this work. The solid lines indicate the best fit relation obtained with Eq. 14. Points and lines are color-coded as a function of time, as indicated in the panel. For clarity, the MS relations are artificially displaced by 0.4 dex from one another. The bottom panel shows the residual distribution as a function of stellar mass. The dashed red line shows the 0 level, corresponding to the best fit value. The red points indicate the mean value of the residuals, while the dotted and dashed-dotted lines indicate the 1 and 3 𝜎 regions, respectively.  \n<!-- image -->  \nevolves consistently with the results of S14. We point out that it evolves roughly as ' 1 , 𝑧 ' 3 . However, this is only an approximation, because the SFR(z) can not be expressed accurately as a power law of ' 1 , 𝑧 ' . The observed evolution is consistent with previous results in the literature, where there is general consensus on the MS normalization evolving as ' 1 , 𝑧 ' 2 GLYPH<149> 8 GLYPH<0> 3 at least up to z GLYPH<24> 4 (e.g. Sargent et al. 2012; Schreiber et al. 2015; Davidzon et al. 2017). Such evolution has been interpreted as reflecting the evolution of the molecular gas mass density and of the consequent availability of gas supply for the galaxy star formation process (Daddi et al. 2008, 2010; Tacconi et al. 2010, 2018).", '3.2 The MS turn-over': "The approach of S14 is very effective in providing a good functional form for the MS. However, it is not trivial to physically interpret Eq.10. For this reason, we explore a different approach by using the fitting function of Lee et al. (2015):  \n𝑆𝐹𝑅 ' 𝑀 ★ ' = 𝑆𝐹𝑅 𝑚𝑎𝑥 GLYPH<157> ' 1 , ' 𝑀 ★ GLYPH<157> 𝑀 0 ' GLYPH<0> 𝛾 ' (11)  \nUnlike polynomial fits, as Eq. 10, the parameters of this model allow to quantify the interesting characteristics of the relation between stellar mass and SFR: i) 𝛾 , the power-law slope at low stellar masses, ii) 𝑀 0 , the turnover mass, and iii) 𝑙𝑜𝑔𝑆𝐹𝑅 𝑚𝑎𝑥 , the maximum value of 𝑙𝑜𝑔𝑆𝐹𝑅 that the function asymptotically approaches at high stellar masses.  \nTo capture the evolution of the MS through Eq. 11, we first fit the MS in bins of time to estimate the time dependence. This allows to obtain 𝑆 0 ' 𝑡 ' , 𝑀 0 ' 𝑡 ' and 𝛾 ' 𝑡 ' to find the best functional form, which includes the time dependence in Eq. 11. The left panel of Fig. 3 shows the results of the best-fit parameters as a function of time.  \nWe find that the exponent 𝛾 does not show any time dependence and it is consistent with the value GLYPH<0> 1. Both 𝑙𝑜𝑔𝑆𝐹𝑅 𝑚𝑎𝑥 and 𝑙𝑜𝑔𝑀 0 depend linearly on time. Thus, they can be expressed as:  \n𝑙𝑜𝑔𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 ' = 𝑎 0 , 𝑎 1 𝑡 (12)  \nand  \n𝑙𝑜𝑔𝑀 0 ' 𝑡 ' = 𝑎 2 , 𝑎 3 𝑡 (13)  \nBy taking into account such evolution, we can write Eq. 11 as a function of time as:  \n𝑙𝑜𝑔𝑆𝐹𝑅 ' 𝑀 ★ GLYPH<150> 𝑡 ' = 𝑎 0 , 𝑎 1 𝑡 GLYPH<0> 𝑙𝑜𝑔 ' 1 , ' 𝑀 ★ GLYPH<157> 10 𝑎 2 , 𝑎 3 𝑡 ' GLYPH<0> 𝑎 4 ' (14)  \nWe use this functional form to fit the MS estimate dataset as a function of mass and time without any binning. The results of the best fitting procedure are indicated in Table 2. The right panel of Fig. 3 shows the MS data points as a function of time with the best fitting curves. The scatter around the best fit is 0.09 dex, indicating that also Eq. 14 provides an excellent fitting function for the evolution of the MS. Consistently with the results obtained with the functional form given by Eq. 10, the MS bends at the high mass end with a turn-over mass that is evolving with time. The exponent 𝛾 does not evolve with time and it is consistent with the value 1. This implies that the MS is well represented by the functional form:  \n𝑆𝐹𝑅 ' 𝑀 ★ ' = 𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 'GLYPH<157> ' 1 , ' 𝑀 0 ' 𝑡 'GLYPH<157> 𝑀 ★ '' (15)  \nwhich is regulated by only two parameters, 𝑀 0 ' 𝑡 ' and 𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 ' . Weestimate the best-fit parameters also for this functional form. The results are included in Table 2. The evolution of the MS shape is regulated by the change of the turn-over mass 𝑀 0 ' 𝑡 ' . This changes of only 25% over the past 9-10 Gyrs, and is a factor of 2 larger in the  \nFigure 4. Comparison of the two MS best fits as a function of time obtained through Eq. 10 (dotted lines) and Eq. 14 (solid lines). All lines are colorcoded as a function of the Universe age as indicated in the figure from 1 to 13 Gyrs.  \n<!-- image -->  \nfirst 3-4 Gyrs. The normalization of the relation, set by 𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 ' parameter, exhibits the strongest evolution and it is consistent within 1 𝜎 with the value found for the normalization in Eq. 10.  \nFig. 4 shows the comparison of the best fit obtained with Eq. 10 and Eq. 14. The two fits exhibit a very consistent evolution in normalization as a function of time. We do observe a discrepancy only at the very low mass and high mass ends, in particular at early epochs, where the uncertainty and scatter of the data-points are the highest. Nevertheless, the two fits are consistent with each other within the 1 𝜎 uncertainty.", '3.3 Towards a physical explanation of the MS shape evolution': "Given Eq. 15, 𝑀 0 ' 𝑡 ' can be interpreted as the mass thresholds between two regimes. At M ★ GLYPH<159> GLYPH<159> 𝑀 0 ' 𝑡 ' the sSFR of galaxies is nearly constant as a function of stellar mass and equal to GLYPH<24> 𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 'GLYPH<157> 𝑀 0 ' 𝑡 ' . At M ★ ¡¡ 𝑀 0 ' 𝑡 ' , the sSFR is progressively suppressed as it is approximated by 𝑆𝐹𝑅 𝑚𝑎𝑥 ' 𝑡 'GLYPH<157> 𝑀 ★ . Thus, the turn-over stellar mass separates a regime of constant star formation rate per unit of stellar mass from a regime of SFR suppression.  \nPopesso et al. (2019b) point out, exploiting the halo mass catalog of Yang et al. (2007), that at z GLYPH<24> 0the region of the MS (within 3 𝜎 from the relation) is completely dominated by central galaxies. Due to the rather tight correlation between the central galaxy stellar mass and the host halo mass (M ℎ ), this implies that the SFG mean host halo mass is increasing along the MS with M ★ . Up to z GLYPH<24> 1.3, we are able to check if this holds by exploiting the cosmic web catalog of Darvish et al. (2017) in the COSMOS field. This is based on the accurate photometric redshifts of the COSMOS15 catalog (Laigle et al. 2016) to identify clusters, groups, and filaments and assign a membership probability to each galaxy up to 𝑧 GLYPH<24> 1 GLYPH<149> 3 and down to stellar masses of 10 9 𝑀 GLYPH<12> . For galaxies with a high probability to be in groups and clusters, the catalog provides also a classification  \nFigure 5. Evolution of the host halo turn-over mass (orange line) as a function of the Universe Age. The shaded orange region indicates the 1 𝜎 uncertainty. This is estimated by combining the error on 𝑀 0 ' 𝑡 ' and the uncertainty in the central galaxy stellar mass and host halo mass correlation of UNIVERSEMACHINE (Behroozi et al. 2019). The solid black line indicates the evolution of the quenching halo mass of Behroozi et al. (2019). The dashed magenta line shows the evolution of the halo mass threshold between cold and hot accretion regime of Dekel & Birnboim (2006).  \n<!-- image -->  \nin central and satellite galaxies, by identifying as central the most massive system. We use the IR selected galaxy catalog of Popesso et al. (2019a) in the COSMOS field, based on the combination of Herschel and Spitzer MIPSdata, to check the central galaxy fraction in the MS region in the redshift window explored by Darvish et al. (2017). The SFR of each galaxy is given by the combination of the UVandIRcontribution (Popesso et al. 2019b). The MS is identified as the region within 3 𝜎 from the MS relation given by Eq. 14 in several redshift bins. We assume 𝜎 = 0 GLYPH<149> 3 dex. At all redshifts up to 𝑧 = 1 GLYPH<149> 3, the MS region turns out to be dominated by central galaxies, which account for 70% of the galaxy population. To check whether the MS region is dominated by central galaxies also at z ¡ 1.3, we use the predictions from the Illustris TNG hydro-dynamical simulation (Pillepich et al. 2018a). In this case, we use as reference the MS determined as in Donnari et al. (2019) on the same data. Also in this case, as expected, central galaxies account for 70-80% of the galaxy population at 𝑧 ¡ 1 GLYPH<149> 3. Thus, is it plausible to assume that the MS is dominated by central galaxies also in the distant Universe.  \nThis aspect is crucial because it allows us to convert the turnover stellar mass of the SFG MS into a turn-over host halo mass , thanks to the correlation between central galaxy stellar mass and host halo mass (Yang et al. 2007; Behroozi et al. 2013; Behroozi et al. 2019). In particular, we use the results of the empirical model UNIVERSEMACHINE of Behroozi et al. (2019) to convert 𝑀 0 ' 𝑡 ' into host halo turn-over mass 𝑀 ℎ 0 ' 𝑡 ' . In Fig. 5 we plot the evolution of 𝑀 ℎ 0 ' 𝑡 ' as a function of time together with the halo mass quenching threshold derived in Behroozi et al. (2019), and the evolution of the transition mass between cold and hot accretion predicted by the theory of mass accretion as in Dekel & Birnboim (2006). The halo mass quenching threshold is defined as the halo mass above which  \nTable 2. The table lists the best fit parameters of Eq.10, Eq. 14 in the first 2 columns. The last column lists the best fit parameters of Eq. 15, which corresponds to Eq. 14 with the slope 𝑎 4 = 1 ( 𝛾 in Eq. 11).  \nthe fraction of quenched galaxies is larger than 50%. The hot/cold transition, instead, is defined as the halo mass at which the cold gas streams coming from the cosmic web filaments are no longer able to penetrate the halo and feed the central galaxy. The curve of Dekel & Birnboim (2006) depends on the halo temperature, but it predicts that at higher redshift ( 𝑧 ¡ 2 GLYPH<149> 5 GLYPH<0> 3) cold gas streams are still able to penetrate massive hot halos. Thus, it represents the threshold between hot and cold gas accretion onto the central galaxy. Behroozi et al. (2019) discuss that the disagreement between the empirical and the theoretical predictions might originate from the non-inclusion of the effects of black hole feedback in the treatment of accretion. This, indeed, might play an important role in affecting the thermodynamical conditions of the circum-galactic-medium (CGM) of the central galaxy, mainly by injecting large quantities of energy into it (Weinberger et al. 2017; Nelson et al. 2019).  \nThe evolution of 𝑀 ℎ 0 with time is remarkably in agreement with the theoretical model of Dekel & Birnboim (2006). The evolution of 𝑀 ℎ 0 ' 𝑡 ' decreases steeply in the first 4-5 Gyrs of the Universe and it reaches a plateau afterward. This would suggest that the turnover mass of the MS might be indicative of the transition between an environment that efficiently sustains the star formation process of the central galaxies, e.g. through cold gas streams, to one which is hostile to the star formation process due to the suppression of the feeding mechanism of the central galaxy (see also Tacchella et al. 2016; Daddi et al. 2022). This suppression is likely maintained in the hot environments, as suggested by nearly all the most sophisticated hydrodynamical simulations, by the interplay between the hot gas in massive halos and central black hole feedback (Voit et al. 2015, 2017; Weinberger et al. 2017; Nelson et al. 2019).", '4 COMPARISON WITH PREVIOUS RESULTS': 'In this Section, we compare our results on the evolution of the MS with previous observational and theoretical results.  \nWefirstcompareourresults to the ones of S14, keeping in mind that the analysis of S14 is limited to the stellar mass range 10 9 GLYPH<149> 7 GLYPH<0> 10 11 𝑀 GLYPH<12> . In Fig. 6 we plot in magenta the homogenized collection of MS relations in seven different redshift bins. The green solid line marks our best-fit estimates at the given redshift, as expressed in Eq. 10, while the green line indicates the best fit of Eq. 14. The MS relation of S14 is perfectly overlapping our relation in their limited stellar mass range. S14 claim that they do not find a bending of the MS at the high mass end. However, we point out that this could be due to the limited stellar mass range considered in the fitting procedure. Indeed, while a single power law might be a good approximation over this range, the extrapolation to larger masses would largely disagree with the data above 10 11 𝑀 GLYPH<12> .  \nThe bending of the MS has been largely discussed in the liter-  \nature and the redshift and stellar mass at which the relation bends vary largely among the different publications included in our collection, from an MS bending only at 𝑧 GLYPH<159> 2 (e.g. Lee et al. 2018), to a correlation that becomes a power law at 𝑧 ¡ 2 (Schreiber et al. 2015; Tomczak et al. 2016), or an MS constant in shape and evolving only in normalization (Whitaker et al. 2012, 2014; Lee et al. 2015). In this work, we find that the MS bends at all redshifts. However, the bending happens above very large stellar masses at early epochs ( 𝑧 ¡ 3). Works limited to lower stellar masses or to relatively small deep fields can not capture this feature. In addition, some of the previous works focus on different aspects, from the analysis of the low mass slope of the relation in a few cases (e.g. Whitaker et al. 2012, 2014), to the shape and normalization at the high mass end in others (e.g. Sherman et al. 2021). We point out, though, that the bending observed in the MS estimated in this work is the result of the combination of all these different estimates, which all tend to be in agreement within a relatively small scatter (0.08 dex), when brought to a common framework. This, perhaps, might suggest that the reported discrepancies are mainly due to possible biases introduced by a limited stellar mass range, the limited volume of the studied deep fields, possible low number statistics, systematics due to SFR indicator, or a different fitting function. Thus, the approach suggested by S14 and implemented here has the potential of overcoming the limitations of the individual analysis and offers a broader and more complete view of the evolution of the MS over a much larger time interval and stellar mass range.  \nFig. 6 shows also the comparison between our results and the predictions of the Illustris TNG300 simulation (Pillepich et al. 2018a). TNG300 is the largest volume simulated in the suite of Illustris TNG hydrodynamical simulations. The choice of TNG300 is driven by the need to sample a sufficiently large volume to capture the rare giant star-forming galaxies at the high mass end of the MS. TheSFRsareaveragedover200Myrandmeasuredwithinaphysical aperture of 2 𝑅 𝑠𝑡𝑎𝑟 , where 𝑅 𝑠𝑡𝑎𝑟 is the stellar half mass radius (see also Donnari et al. 2019). The shaded region in Fig. 6 shows the distribution of the IllustrisTNG galaxies, color-coded according to the galaxy number density in bins of SFR and 𝑀 ★ . The IllustrisTNG MS is estimated as a running mean as a function of 𝑀 ★ and it is indicated by the orange-yellow line. We also show in yellow the MS of Donnari et al. (2019), which is estimated up to 𝑧 GLYPH<24> 2 by mocking the UVJ selection of Whitaker et al. (2012). The two estimates agree remarkably well up to stellar masses of GLYPH<24> 10 10 GLYPH<149> 5 𝑀 GLYPH<12> and diverge at larger stellar masses, where the UVJ selection of Donnari et al. (2019) exclude most of the low SFR systems.  \nIt is interesting to notice that a bending of the predicted MS is observed also in Illustris TNG at least up to z GLYPH<24> 3. This was already pointed out by Donnari et al. (2019) at 𝑧 GLYPH<159> 2. At higher redshift, the relation is steeper than the observed MS, but this could be due to a  \nFigure 6. MS of star-forming galaxies in several redshift bins from 𝑧 GLYPH<24> 6 to 𝑧 GLYPH<24> 0. The shaded blue region in each panel indicates the distribution of simulated Illustris TNG galaxies in the SFR-stellar mass plane. The blue scale is according to the number density of galaxies. The orange line shows the running median in bins of stellar mass for the simulated galaxies, while the connected yellow triangles show the relation of Donnari et al. (2019) limited to the UVJ selected galaxies in Illustris TNG. The magenta points show the data collected in this paper. The green and purple solid lines indicate our best fit according to Eq. 10 and 14.  \n<!-- image -->  \nFigure 7. Evolution of sSFR as a function of redshift. The black lines show the sSFR evolution given by the best fit of Eq. 14 at three stellar mass values: 10 9 (dot-dashed line), 10 10 (solid line) and 10 11 𝑀 GLYPH<12> (dashed line). We provide the same relation for the S14 best fit at 10 10 (red line). We plot for comparison the data of several works in literature with symbols indicated in the figure. The yellow, green, purple and light purple solid lines show the predictions of EAGLE simulation of Furlong et al. (2015), The MUNICH simulation of Henriques et al. (2015), the Illustris TNG as estimated here and the one of Donnari et al. (2019), respectively. The magenta solid line shows the best fitting model of Behroozi et al. (2019).  \n<!-- image -->  \nlack of massive star-forming galaxies in the simulations, as pointed out below.  \nWe notice a good agreement between the IllustrisTNG predictions and our results at early epochs ( 𝑧 GLYPH<24> 6). The observed MS lies over the simulated relation. However, no giant star-forming galaxies are observed at all in the TNG300 volume at 𝑧 GLYPH<24> 6 above 10 10 GLYPH<149> 5 𝑀 GLYPH<12> . The observations, instead, indicate that such galaxies populate the high mass end of the relation, in agreement with the results of the GSMF of active galaxies of Davidzon et al. (2017). This is also confirmed by the recent JWST discoveries of very massive blue galaxies in the distant Universe (e.g. Santini et al. 2022; Finkelstein &Bagley 2022; Castellano et al. 2022; Donnan et al. 2022)  \nMassive star-forming galaxies appear in the TNG simulation only between 𝑧 GLYPH<24> 4 to GLYPH<24> 3, with a larger SFR compared to observations, in particular at 𝑧 GLYPH<24> 4. At this epoch the predicted and observed low mass slope (below stellar masses of 10 10 GLYPH<149> 5 𝑀 GLYPH<12> ) are in agreement within 1 𝜎 . At later epochs ( 𝑧 GLYPH<159> 2), the predicted MS is systematically below the observations by GLYPH<24> 0 GLYPH<149> 5 dex, according to the results of Donnari et al. (2019). Observations and predictions are again in agreement at 𝑧 GLYPH<24> 0 up to GLYPH<24> 5 GLYPH<2> 10 10 𝑀 GLYPH<12> . Above this mass, no more SFGs are detected in the simulations at odds with observations.  \nThe tension between observations and simulations in the evolution of the MS is not new and has been extensively discussed in the literature (see, for instance, Katsianis et al. 2020; Nelson et al. 2021). From the observational point of view, the use of sophisticated SED fitting codes, such as the Prospector or its incarnation of Prospector𝛼 used in Leja et al. (2022), might apparently solve this  \nissue. Indeed, adding an extra population of old stars might justify an overestimation of the previous Spitzer and Herschel based SFR estimates and reconcile ad hoc the observed and simulated MS in specific redshift bins. However, it is worth pointing out that such an approach might simply move the tension somewhere else. Indeed, the stellar masses provided by Prospector increase considerably with respect to previous works. This increase would shift the stellar mass function leading to a substantial disagreement with the predicted ones. These, indeed, are overall in agreement with previous measurements, as shown, for instance, in Pillepich et al. (2018b) for Illustris TNG and in the detailed comparison of Thorne et al. (2021).  \nIt is, perhaps, more interesting to point out that the largest disagreement between the observed and the simulated MS is in the distribution of the star-forming galaxy population along the MS. Fig. 6 shows quite clearly that, in the observations, star-forming galaxies with stellar masses above GLYPH<24> 5 GLYPH<2> 10 10 𝑀 GLYPH<12> exhibit a much slower evolution than predicted by the simulations. Indeed, they should have formed by 𝑧 GLYPH<24> 4, while in Illustris TNG they do not appear at all. Furthermore, the high mass end of the MS, above this mass threshold, is almost completely evacuated by 𝑧 GLYPH<24> 0 in the simulation. Instead, many works in the literature, based on the SDSS galaxy spectroscopic sample or the WISE survey, show clearly that this region of the SFR-stellar mass plane is highly populated up to stellar masses of 3 GLYPH<2> 10 11 𝑀 GLYPH<12> (see, for instance, Popesso et al. 2019b, for a collection of the available datasets in the local Universe). Such discrepancy can not be ascribed to an overestimation of the observed SFR, as recently proposed (Leja et al. 2022). Indeed, this would simply lead to a discrepancy of the MS normalization rather than to the lack of the massive SFG sub-population. Thus, it is, perhaps, more likely that the predicted fast evolution of massive SFG is due to an over-efficient SF quenching process in simulations. We point out that the use of the larger stellar masses estimated by Prospector would increase the observed discrepancy at the high mass end, not only in the local Universe but also at higher redshift (0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 2).  \nWe provide in Fig. 7 a comparison of the sSFR evolution with simulations, as well as other works in literature. We plot the evolution of the sSFR derived from the best fit of Eq. 14 computed in three different stellar mass bins: 10 9 (black dot-dashed line), 10 10 (black solid line) and 10 11 𝑀 GLYPH<12> (black dashed line). The red line shows the result of S14 at 10 10 𝑀 GLYPH<12> , which is in agreement with our result at the same stellar mass. The prediction of IllustrisTNG, up to 𝑧 GLYPH<24> 6, those from Donnari et al. (2019), up to 𝑧 GLYPH<24> 2, as well as the best fit of the empirical model UNIVERSEMACHINE (Behroozi et al. 2019), all computed for M ★ = 10 10 𝑀 GLYPH<12> , are shown in purple, lavender and magenta, respectively. The simulated data of Illustris TNG are in agreement with observations at high redshift and in the local Universe, but the slope of the predicted sSFR evolution is steeper than observations, as pointed out in the previous paragraph. The redshift evolution of sSFR obtained from Behroozi et al. (2019) is perfectly overlapping with observations at the same stellar mass. Finally, we show in Fig. 7 the redshift evolution of the sSFR obtained from the EAGLE hydrodynamical simulation (Furlong et al. 2015) and the Munich simulation of Henriques et al. (2015). Both predictions lie below the observations and exhibit a steeper relation. However, we point out that, as discussed by Davidzon et al. (2018), such comparison is complicated by the fact that simulated galaxies are not selected to be SFGs. Thus, the slope and normalization of the sSFR-redshift relation is biased by quiescent galaxies, which might have a faster evolution.  \nWeconcludethat the star formation activity of simulated galaxies is declining much faster than in observations between redshift  \nGLYPH<24> 3 and GLYPH<24> 0. This leads to an MS normalization too low with respect to observations in the same redshift window. The disagreement is more significant for the most massive systems above 10 11 𝑀 GLYPH<12> , which completely disappear from the simulated MS in the local Universe, at odds with observations. We do not find evidence for a clear overestimation or underestimation of the SFRs in specific redshift bins, that might solve the tension. As discussed in Donnari et al. (2019), the observed discrepancy might point to some fundamental limitations in our understanding of the processes governing the star formation evolution in galaxies, such as the role of cold and hot accretion and of supernovae and black hole feedback.', '5 SUMMARY AND CONCLUSION': 'We compile a collection of the most important publications regarding the evolution of the MS of star forming galaxies in the widest range of redshift (0 GLYPH<159> 𝑧 GLYPH<159> 6), stellar mass (10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> ) and SFR (0 GLYPH<0> 500 𝑀 GLYPH<12> 𝑦𝑟 GLYPH<0> 1 ) ever probed. We convert all observations to a common calibration to check for consistency in the literature estimates and to study the evolution of the relation with different approaches. We find a remarkably good agreement between the different estimates at any stellar mass and time. The resulting MS exhibits a curvature towards the high stellar masses, which is slowly evolving with time. We provide two functional forms, which take into account the time evolution of the MS normalization and slope. Following the approach of S14, we estimate the best polynomial fitting form in the logSFR GLYPH<0> logM ★ space, by studying the evolution with time of the logSFR at fixed stellar mass. A second-order polynomial form is well representing the relation, with a non-evolving quadratic term. The normalization is evolving as a power law of the Universe age. The slow evolution of the linear term reproduces the steepening of the relation towards the first 3-4 Gyrs of the Universe. We provide, as alternative fitting form, the one of Lee et al. (2015), which has the advantage of being expressed by physical parameters. These are the low mass slope, the normalization, and the turn-over mass ( 𝑀 0 \' 𝑡 \' ). While the slope does not evolve with time, normalization and turn-over mass evolve as a power law of the Universe age. The turn-over mass, in particular, determines the MS shape. It marginally evolves with time, and it is responsible for the steepening of the relation towards 𝑧 GLYPH<24> 4 GLYPH<0> 6. At stellar masses below 𝑀 0 \' 𝑡 \' , SFGs have a constant sSFR, while above 𝑀 0 \' 𝑡 \' the sSFR is suppressed. As the MS region is dominated by central galaxies, we use the relation between central galaxy stellar mass and host halo mass to convert 𝑀 0 \' 𝑡 \' into a "turn-over host halo mass". We find that its evolution is remarkably consistent with the one of the halo mass threshold between cold and hot accretion regimes predicted by Dekel & Birnboim (2006). This might indicate that 𝑀 0 \' 𝑡 \' defines the transition between an environment able to sustain the SF process of the central galaxy, for instance through cold gas streams as predicted by Dekel & Birnboim (2006), to a regime hostile to the same process. The latter might be the result of the interplay between the hot gas in massive halos and the black hole feedback generated by the central galaxy itself.  \nThe comparison of our results with the state-of-the-art hydrodynamical simulations shows that the simulated MS shape is qualitatively consistent with the observations. However, the normalization of the simulated relation is systematically lower by at least 0.2 to 0.5 dex with respect to observations at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3. As a consequence, the sSFR of SFGs evolves more rapidly in simulated galaxies with respect to the observed ones. We do not find clear evidence for an overestimation of the observed SFR at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3  \nwith respect to simulated SFR. This might suggest that the feedback implemented in the simulations is likely too efficient in suppressing the star formation activity of galactic systems, in particular, above 5 GLYPH<2> 10 10 𝑀 GLYPH<12> (see also Katsianis et al. 2021; Corcho-Caballero et al. 2021).', 'ACKNOWLEDGEMENTS': 'We would like to deeply thank the referee, A. Katsianis, for helping in taking a different perspective and for the enormous help in catching up with the literature. This paper was ready for submission at the end of 2019, but P.P. had to stop working for about 2 years because of the pandemic. The referee helped significantly in improving the manuscript and we consider this as one of the most positive and constructive experiences of a review process.', 'DATA AVAILABILITY': 'The data underlying this article are publicly available in the articles listed in Table 1 and in their online supplementary material.', 'APPENDIX A: THE SELECTED MS': 'Here we list and describe the MS estimates included in the analysis and the data they are based on:  \n- GLYPH<15> Speagle et al. (2014) do not provide the collection of data used for calibrating the average MS estimates but only the best fits. A large number of best fits are provided based on different subsamples, depending on the star-forming galaxy selection method, and on the range of stellar masses and cosmic time included in the fit. We use in this paper the fit n. 64 based on a "mixed selection", which S14 consider as more inclusive of the SFG population (see Section 2.3 for a detailed discussion). The fit is restricted to stellar masses between 10 9 GLYPH<149> 7 𝑀 GLYPH<12> and 10 11 𝑀 GLYPH<12> and it excludes the first 2 Gyrs of the Universe Age. In order to properly populate the S14 MS given by the fit n. 64, we use the data of Fig n. 4, which shows the value of the SFR as a function of time in 4 stellar mass bins. We distribute randomly the data within the stellar mass bins and use the time information to estimate the SFR as a function of stellar mass and time in the range 10 9 GLYPH<149> 7 GLYPH<0> 10 11 𝑀 GLYPH<12> and 2-12 Gyrs of Universe age. The SFRs obtained in this way are calibrated to the KE12 calibration and a Kroupa IMF.\n- GLYPH<15> Rodighiero et al. (2014) provide the MS estimate of BzK selected galaxies in the redshift bin 1 GLYPH<149> 4 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 5 in the COSMOS field, on the basis of various SFR indicators including UV emission, 𝐻𝛼 emission, MIR and FIR emission. In this analysis, we use the four data points obtained through the stacking analysis of the BzK sample in the COSMOS PACS maps in the stellar mass range 10 10 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> . The SFR and stellar masses are obtained with a Salpeter IMF and are corrected to the KE12 calibration and a Kroupa IMF.\n- GLYPH<15> Heinis et al. (2014) provide the MS based on the stacking UV selected SFGs in bins of FUV luminosity in the HerMES maps (Oliver et al. 2012). The derived SFR is based on the combination of NUV and IR luminosities, with a Chabrier IMF. The MS is  \nestimated between 10 9 GLYPH<149> 5 and 10 11 GLYPH<149> 3 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 1 GLYPH<149> 5, 𝑧 GLYPH<24> 3 and 𝑧 GLYPH<24> 4. We include in our analysis the stacked points at the observed stellar masses and redshift given in the paper after correcting them to the KE12 calibration.  \n- GLYPH<15> Whitaker et al. (2014) provide the MS based on stacking of UVJ selected star-forming galaxies on Spitzer MIPS 24 𝜇 mdata in the 3D-HST CANDELS fields. The derived SFR is based on the combination of NUV and IR luminosities, with a Chabrier IMF. The considered redshift range is 0.5-2.5 and the MS is estimated in the following stellar mass and redshift ranges: 10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 2 𝑀 GLYPH<12> at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 1, 10 9 GLYPH<0> 10 11 GLYPH<149> 3 𝑀 GLYPH<12> at 1 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 5, 10 9 GLYPH<149> 2 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> at 1 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 2 and 10 9 GLYPH<149> 3 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> at 2 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 5. We include in our analysis the stacked points at the observed stellar masses and redshift given in the paper, recalbrated to KE12.\n- GLYPH<15> Chang et al. (2015) provide SFR based on SED fitting results of MAGHPHYS code, from GALEX to WISE data in the local Universe. As already reported by Popesso et al. (2019b), the catalog provides SFRs underestimated with respect to the H 𝛼 based and IR-based SFR estimates. Salim et al. (2016) discuss that this is likely due to the fact that the MAGHPHYS SED fitting results are mostly driven by the higher SNR 12 𝜇 m WISE data-point than the low SNR 22 𝜇 mWISE data, leading to artificially low SFR. Chang et al. (2015) report that their MS at 𝑧 GLYPH<159> 0 GLYPH<149> 1 is systematically lower (0.15 dex) with respect to other H 𝛼 and IR based MS. We estimate the MS by using the best fit relation in stellar mass bins of 0.15 dex to take into account stellar mass uncertainties and after correcting the SFRs for the 0.15 dex offset, which brings them in agreement to the Ke12 calibration.\n- GLYPH<15> Lee et al. (2015) use a ladder of SFR indicators to study the MS in the COSMOS field at 0 GLYPH<149> 3 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 3. Due to the flux limit of the Spitzer MIPS 24 𝜇 m and Herschel PACS and SPIRE catalogs, the combination of NUV and IR luminosities provides the SFR only for highly star-forming objects (starburst and massive SFGs). For less active or less dusty objects they use dust-corrected NUV luminosities (see Fig. 3 of the cited paper for the contribution of any SFR indicator as a function of stellar mass). Star-forming galaxies are selected according to their (M 𝑁𝑈𝑉 -M 𝑅 ) GLYPH<0> (M 𝑅 -M 𝐽 ) colors with a Chabrier IMF. The MS is estimated as the median SFR value in equally populated bins of stellar mass in the following ranges: 10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 𝑀 GLYPH<12> at 𝑧 = 0 GLYPH<149> 36, 10 9 GLYPH<0> 10 11 GLYPH<149> 𝑀 GLYPH<12> at 𝑧 = 0 GLYPH<149> 55, 10 9 GLYPH<0> 10 11 GLYPH<149> 𝑀 GLYPH<12> at 𝑧 = 0 GLYPH<149> 70, and between 10 9 GLYPH<149> 3 and 10 11 GLYPH<149> 2 𝑀 GLYPH<12> at 𝑧 = 0 GLYPH<149> 85, 𝑧 = 0 GLYPH<149> 99 and 𝑧 = 1 GLYPH<149> 19. We include in our analysis the median points at the observed stellar masses and redshift, after correcting the median into the mean value and to the KE12 calibration.\n- GLYPH<15> Ilbert et al. (2015) use the flux-limited sample of Spitzer MIPS 24 𝜇 mand Herschel PACSandSPIREoftheCOSMOSfieldtostudy the distribution of galaxies in the SFR-stellar mass plane. Galaxies are further selected according to their (M 𝑁𝑈𝑉 -M 𝑅 ) GLYPH<0> (M 𝑅 -M 𝐽 ) colors. The SFR is estimated with a combination of NUV and IR luminosities. The MS is identified as the peak of the SFR distribution at several stellar mass bins at stellar masses above 10 10 𝑀 GLYPH<12> and in the redshift range 0 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 4. This corresponds to the median SFR in a log-normal distribution. We include in our analysis the median points at the observed stellar masses and redshift given in the paper, after correcting the median to the mean and to the KE12 calibration.\n- GLYPH<15> Tasca et al. (2015) use the spectroscopically selected sample of the VIMOS Ultra Deep Survey (VUDS, Le Fèvre et al. 2015). Galaxies are selected from a combination of photometric redshifts, as well as from color selection criteria like LBG, combined with a fluxlimit 22 GLYPH<149> 5 GLYPH<159> 𝑖 𝐴𝐵 GLYPH<159> 25. The SFR is derived via SED fitting techniques with Le Phare , assuming the BC03 SPS model, a Chabrier  \nIMF and a Calzetti et al. (2000) extinction law. The MS is provided in the following stellar mass and redshift ranges: 10 7 GLYPH<149> 5 GLYPH<0> 10 10 𝑀 GLYPH<12> at 0 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 7, 10 8 GLYPH<149> 5 GLYPH<0> 10 11 𝑀 GLYPH<12> at 0 GLYPH<149> 7 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 5, 10 9 GLYPH<0> 10 11 𝑀 GLYPH<12> at 1 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 5 and at 2 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 5, 10 9 GLYPH<149> 3 GLYPH<0> 10 10 GLYPH<149> 6 𝑀 GLYPH<12> at 3 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 4 GLYPH<149> 5 and 10 9 GLYPH<0> 10 11 𝑀 GLYPH<12> at 𝑧 ¡ 4 GLYPH<149> 5. We include in our analysis the mean points at the observed stellar masses and redshift, corrected to the KE12 calibration.  \n- GLYPH<15> Salmon et al. (2015) study the evolution of the slope and scatter of the SFR-stellar mass relation for galaxies at 3 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 6 GLYPH<149> 5 using multi-wavelength photometry in GOODS-N from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) and Spitzer Extended Deep Survey. They use an updated, Bayesian spectral-energy distribution fitting method that incorporates effects of nebular line emission, star formation histories that are constant or rising with time, and different dust attenuation prescriptions (starburst and Small Magellanic Cloud). They use a modified version of Bruzual & Charlot (2003) stellar population synthesis models, and a Salpeter IMF.\n- GLYPH<15> Renzini & Peng (2015) provide a fit of the local MS at 𝑧 GLYPH<159> 0 GLYPH<149> 085 in the stellar mass range 10 9 GLYPH<0> 10 10 GLYPH<149> 5 𝑀 GLYPH<12> . The SFR is based on dust-corrected SDSS H 𝛼 fluxed for systems classified as star forming in the BPT diagram and through the D4000 break for non-active systems or AGN hosts (Brinchmann et al. 2004). SFR and stellar masses are estimated with a Chabrier IMF. The MS is identified as the peak of the SFR distribution at fixed stellar mass in the SFR-stellar mass plane. This corresponds to the median SFR in a log-normal distribution. We estimate the MS by using the best fit relation in stellar mass bins of 0.15 dex to take into account stellar mass uncertainties. The MS estimates are included after correcting the median into the mean value and to the KE12 calibration.\n- GLYPH<15> Schreiber et al. (2015) perform a stacking analysis of UVJ selected galaxies in the deep Herschel PACSmapsoftheCANDELS fields. The mean IR luminosity derived for the stacks is combined with the mean FUV luminosity to derive the SFR. SFR and stellar masses are estimated with a Salpeter IMF. The MS is given in the following stellar mass and redshift ranges: 10 9 GLYPH<0> 10 11 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 0 GLYPH<149> 5, 𝑧 GLYPH<24> 1 and 𝑧 GLYPH<24> 1 GLYPH<149> 6, 10 9 GLYPH<149> 8 GLYPH<0> 10 11 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 2 and 𝑧 GLYPH<24> 3. We include in our analysis the stacked points at the observed stellar masses and redshift, corrected to the KE12 calibration.\n- GLYPH<15> de los Reyes et al. (2015) use GLYPH<24> 300 H 𝛼 -selected galaxies at 𝑧 GLYPH<24> 0 GLYPH<149> 8 to study the MS. They use deep optical spectra obtained with the IMACS spectrograph at the Magellan telescope to measure strong oxygen lines. They combine spectral information with rest-frame UV-to-optical imaging, which allows them to determine stellar masses and dust attenuation corrections, and H 𝛼 narrow-band imaging, which provides a robust measurement of the instantaneous SFR. The SEDs are fit with a library of Bruzual & Charlot (2003) stellar population synthesis models. The model libraries are built with a wide range of SFHs and metallicities, as described in Salim et al. (2007), and updated da Cunha et al. (2008). Each model is attenuated according to the prescription of Charlot & Fall (2000). The dust attenuation in the SED fitting is mainly constrained by the UV slope.  \nTheir sample spans stellar masses of GLYPH<24> 10 9 ˘6 GLYPH<2> 10 11 𝑀 GLYPH<12> . They find consistency with previous results about the MS at similar redshift. We do not implement any correction in stellar mass and SFR.  \n- GLYPH<15> Erfanianfar et al. (2016) provide the MS at 𝑧 GLYPH<159> 1 GLYPH<149> 1. The SFR is based on the flux-limited sample of Spitzer MIPS 24 𝜇 m and Herschel PACS and SPIRE available in the ECDFS and COSMOS fields. The SFR is derived from the far-IR flux with a Chabrier IMF. The MS is retrieved via 𝜎 -clipping. It is estimated in the stellar mass range 10 9 GLYPH<0> 10 11 GLYPH<149> 5 at 0 GLYPH<149> 15 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 5 and 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 1.', '16': "We include in our analysis the stacked points at the observed stellar masses and redshift, corrected to the KE12 calibration.  \n- GLYPH<15> Tomczak et al. (2016) perform a stacking analysis of UVJ selected galaxies in the deep Herschel maps of the CANDELS fields. The mean IR luminosity derived for the stacks is combined with the mean NUV luminosity to derive the SFR. SFR and stellar masses are estimated with a Chabrier IMF. The MS is given in the following stellar mass and redshift ranges: 10 8 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 2 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 0 GLYPH<149> 6, 𝑧 GLYPH<24> 0 GLYPH<149> 9, 𝑧 GLYPH<24> 1 GLYPH<149> 1 and 𝑧 GLYPH<24> 1 GLYPH<149> 4, 10 8 GLYPH<149> 7 GLYPH<0> 10 11 GLYPH<149> 2 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 1 GLYPH<149> 7, 10 9 GLYPH<0> 10 11 GLYPH<149> 2 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 2 GLYPH<149> 25 and 𝑧 GLYPH<24> 2 GLYPH<149> 75, 10 9 GLYPH<149> 5 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> at 𝑧 GLYPH<24> 3 GLYPH<149> 5. We include in our analysis the stacked points at the observed stellar masses and redshift, corrected to the KE12 calibration.\n- GLYPH<15> Santini et al. (2017) derive the MS in the ultra-deep Hubble Space Telescope Frontier fields. The MS is derived via SED fitting with a Salpeter IMF, BC03 SPS model and a Calzetti et al. (2000) extinction law. The MS location is derived via 𝜎 -clipping at the following stellar masses and redshift ranges: 10 8 GLYPH<0> 10 10 GLYPH<149> 6 𝑀 GLYPH<12> at 1 GLYPH<149> 3 GLYPH<159> 𝑧 GLYPH<159> 2, 10 8 GLYPH<0> 10 11 𝑀 GLYPH<12> at 2 GLYPH<159> 𝑧 GLYPH<159> 3 and 3 GLYPH<159> 𝑧 GLYPH<159> 4, 10 8 GLYPH<0> 10 11 at 4 GLYPH<159> 𝑧 GLYPH<159> 5, and 10 8 GLYPH<149> 6 GLYPH<0> 10 11 𝑀 GLYPH<12> at 5 GLYPH<159> 𝑧 GLYPH<159> 6. We estimate the MS by using the best fit relations in stellar mass bins of 0.25 dex to take into account stellar mass uncertainties. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Kurczynski et al. (2016) utilize photometry in the Hubble Ultradeep Field (HUDF12) and Ultraviolet Ultra Deep Field (UVUDF) campaigns and CANDELS/GOODS-S to estimate the SFR via SED fitting. They use a Salpeter IMF, BC03 SPS model and a Calzetti et al. (2000) extinction law. The MS is identified via 𝜎 -clipping between 10 7 and 10 11 𝑀 GLYPH<12> in the following redshift ranges: 0 GLYPH<149> 5 GLYPH<0> 1 GLYPH<149> 0, 1 GLYPH<149> 0 GLYPH<0> 1 GLYPH<149> 5, 1 GLYPH<149> 5 GLYPH<0> 2 GLYPH<149> 0, 2 GLYPH<149> 0 GLYPH<0> 2 GLYPH<149> 5 and 2 GLYPH<149> 5 GLYPH<0> 3 GLYPH<149> 0. We point out that beyond 𝑧 GLYPH<24> 1 GLYPH<149> 5, the stellar mass range 10 10 GLYPH<149> 2 GLYPH<0> 10 11 is populated by less than 20-25 systems. Thus, we limit the use of the MS estimates to the 10 7 GLYPH<0> 10 10 GLYPH<149> 2 𝑀 GLYPH<12> stellar mass range. We estimate the MS by using the best fit relations in stellar mass bins of 0.25 dex to take into account stellar mass uncertainties. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Pearson et al. (2018) use the SED modeling and fitting tool CIGALE to generate flux density priors in the Herschel SPIRE bands in the COSMOS field. These priors are fed into a deblending tool, called 𝑋𝐼𝐷 , , to extract flux densities from the SPIRE maps. As the last step, multi-wavelength data are combined with the extracted SPIRE flux densities to constrain SEDs and provide stellar masses and SFRs. These are used to populate the SFR𝑀 GLYPH<3> plane over the redshift range 0 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 6. The SED fitting is performed with a Chabrier IMF, BC03 SPS model, and a Charlot & Fall (2000) extinction law. Star-forming galaxies are selected through the UVJ selection as in Whitaker et al. (2014). The MS is estimated in the following stellar mass and redshift ranges: between 10 9 and 10 11 𝑀 GLYPH<12> at 0 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 5 and 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 8, between 10 9 GLYPH<149> 5 and 10 11 𝑀 GLYPH<12> at 0 GLYPH<149> 8 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 1, 1 GLYPH<149> 1 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 4 and 1 GLYPH<149> 4 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 8, between 10 10 and 10 11 𝑀 GLYPH<12> at 1 GLYPH<149> 8 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 3 and 2 GLYPH<149> 3 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 9, between 10 10 GLYPH<150> 5 and 10 11 𝑀 GLYPH<12> at 2 GLYPH<149> 9 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 9 and 3 GLYPH<149> 9 GLYPH<159> 𝑧 GLYPH<159> 4 GLYPH<149> 9 and above 10 11 𝑀 GLYPH<12> at 4 GLYPH<149> 9 GLYPH<159> 𝑧 GLYPH<159> 6. We estimate the MS by using the best fit relations in stellar mass bins of 0.25 dex to take into account stellar mass uncertainties. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Belfiore et al. (2018) provide a measure of the local MS at 𝑧 GLYPH<159> 0 GLYPH<149> 05 in the MaNGA sample (Bundy et al. 2015). Star-forming galaxies are selected according to the BPT diagram in any spaxel of the MaNGA map. The SFR is estimated through the integrated, extinction-corrected H 𝛼 emission over the galaxy region, converted via Robert C Kennicutt (1998) law and with a Chabrier IMF. The MS is estimated as the mean SFR in bins of stellar masses in the range 10 9 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> . We include in our analysis the mean points  \nat the observed stellar masses and redshift and correct to the Kroupa IMF.  \n- GLYPH<15> Davidzon et al. (2018) estimate the evolution of the specific star formation rate at a fixed mass of 10 10 GLYPH<149> 3 𝑀 GLYPH<12> . The SFR is derived by the evolution of the galaxy stellar mass function of active galaxies of Davidzon et al. (2017), under the assumption that the stellar mass growth of a galaxy is mainly driven by the integral of the SFR history, minus the quantity of stellar mass returned to the interstellar medium. The GSMF of active galaxies, and thus, the derived SFR evolution is based on the (M 𝑁𝑈𝑉 -M 𝑅 ) GLYPH<0> (M 𝑅 -M 𝐽 ) colors. The SED fitting technique used in this analysis is based on a Chabrier IMF. The MS is estimated at 𝑧 GLYPH<24> 2.22, 2.75, 2.25, 3.75, 5, 6 and 7. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Lee et al. (2018) use the deep CANDELS observations in the GOODS North and South fields to revisit the correlations between stellar mass and star formation rate in galaxies at 1 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 4. The quantities are estimated via SED fitting with a Chabrier IMF and a BC03 SPS model. SFGs are selected through the UVJ selection. The MS is estimated at stellar masses above 10 9 𝑀 GLYPH<12> and in the following redshift bins: 1 GLYPH<149> 2 GLYPH<0> 1 GLYPH<149> 5, 1 GLYPH<149> 5 GLYPH<0> 2, 2 GLYPH<0> 2 GLYPH<149> 8 and 2 GLYPH<149> 8 GLYPH<0> 4. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Iyer et al. (2018) use star formation histories (SFHs) reconstructed via the Dense Basis method of Iyer & Gawiser (2017) at 0.5 < z < 6 in the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey GOODS-S field to study the nature and evolution of MS. They use the reconstructed SFHs as trajectories in SFR-stellar mass plane. This allows them to study galaxies at epochs earlier than observed by propagating them backward in time along these trajectories. To generate spectra corresponding to a galaxy with a given basis SFH, they use the Flexible Stellar Population Synthesis (FSPS) mode (Conroy et al. 2009; Conroy & Gunn 2010). They use a Chabrier IMF, a Calzetti et al. (2000) attenuation law, and IGM absorption according to Madau et al. (1996) prescription. They study the MS at z = 1, 2, 3, 4, 5, 6 using both direct fits to galaxies observed at those epochs and 𝑆𝐹𝑅 GLYPH<0> 𝑀 GLYPH<3> trajectories of galaxies observed at lower redshifts. In order to exclude quiescent galaxies, they impose a cutoff, excluding systems that are at a distance larger than 0.4 dex from the best fit MS correlation. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Popesso et al. (2019b) provide measures of the local MS at 𝑧 GLYPH<159> 0 GLYPH<149> 085 and stellar masses ¡ 10 10 𝑀 GLYPH<12> based on several SFR indicators. These include SDSS dust corrected H 𝛼 emission (Brinchmann et al. 2004), WISE 22 𝜇 memission, dust corrected UV emission (Salim et al. 2016) and Herschel PACS and SPIRE emission at 100, 160, 250, 350 and 500 𝜇 m (Valiante et al. 2016). We include all of them in the analysis. The SFR are all based on a Chabrier IMF. The MS is estimated as the peak of the SFR distribution in several stellar mass bins in the SFR-stellar mass plane. We include in our analysis the mean points at the observed stellar masses and redshift. We correct the MS to the Kroupa IMF.\n- GLYPH<15> Popesso et al. (2019a) measure the evolution of the MS up to z GLYPH<24> 2.5 in the stellar mass range 10 10 GLYPH<0> 10 11 GLYPH<149> 5 𝑀 GLYPH<12> . The SFR is derived by a combination of IR fluxes based on Spitzer MIPS24 𝜇 m and Herschel PACS and SPIRE data and NUV emission available in the CANDELS and COSMOS fields. The SFR and stellar masses are based on a Chabrier IMF. The MS is estimated as the peak of the SFR distribution in several stellar mass and redshift bins at 0 GLYPH<149> 3 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 5, 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 0 GLYPH<149> 8, 0 GLYPH<149> 8 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 2, 1 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 6, 1 GLYPH<149> 6 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 2 and 2 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 2 GLYPH<149> 5. This corresponds to the median SFR in a log-normal distribution. We correct the median into the mean and\n- GLYPH<15> Barro et al. (2019) use a WFC3 F160W (H-band) selected cata-  \nlog in the CANDELS/GOODS-N field containing photometry from the UV to far-infrared, photometric redshifts and stellar parameters derived from the analysis of the multi-wavelength data. Stellar masses are estimated with FAST and are consistent with those of the 3D-HST survey (e.g. Whitaker et al. 2014). They use a ladder of SFR indicators. The SFR ladder consists of three steps that differ on the amount of SFR indicators that are available for each galaxy, namely, UV, mid-IR, and far-IR. The first step is to compute the SFR from the UV luminosity. In this respect, they estimated the total UV luminosity from SED fitting, and use the relation of Robert C Kennicutt (1998), applying a correction for dust attenuation. The UV attenuation is inferred directly from the best-fit model to the overall SED, which assumes a Calzetti et al. (2000) attenuation law. For the second and third steps of the ladder, they obtain the total SFR by adding the contributions from the IR emission, either from mid-IR or from far-IR, when available, and the observed, unobscured UV component. They use the calibration of Robert C Kennicutt (1998) corrected to a Chabrier IMF. The authors state that stellar masses and star formation rates are consistent with previous results. The MS is estimated in the redshift bins 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 0, 1 GLYPH<149> 0 GLYPH<159> 𝑧 GLYPH<159> 1 GLYPH<149> 8, and 1 GLYPH<149> 8 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 0 over the stellar mass range GLYPH<24> 10 8 GLYPH<149> 6 GLYPH<0> 10 11 GLYPH<149> 3 𝑀 GLYPH<12> . It is measured as the running median of the star-forming galaxies identified with the UVJ criterion. We correct the median to the mean and to the Kroupa IMF.  \n- GLYPH<15> Leslie et al. (2020) measure the evolution of the MS in the COSMOS field by exploiting the VLA-COSMOS 3 GHz Large Project dataset. The MS is estimated in the range 0 GLYPH<149> 3 GLYPH<159> 𝑧 GLYPH<159> 6 in the stellar mass range 10 9 GLYPH<0> 10 11 𝑀 GLYPH<12> with different completeness limits as a function of the redshift. The median SFR is estimated via stacking in radio data in several redshift and stellar mass bins. We consider here only the bins where the prior galaxy sample is complete in stellar mass, as indicated by the authors. The SFR is estimated with a Chabrier IMF. SFGs are selected as in Ilbert et al. (2013) through the NUVRJ color-color selection. We correct the median to the mean and to the Kroupa IMF.\n- GLYPH<15> Thorne et al. (2021) apply ProSpect (Robotham et al. 2020) in a parametric mode to multiwavelength photometry from the Deep Extragalactic VIsible Legacy Survey (DEVILS, Davies et al. 2018) in order to measure stellar and dust masses and SFR for galaxies in the COSMOSfield in the 0 GLYPH<159> 𝑧 GLYPH<159> 9 redshift range. They use the new DEVILS 𝑈𝑉 to 𝐹𝐼𝑅 photometry derived using profound (Robotham et al. 2018). ProSpect uses the Bruzual & Charlot (2003) stellar libraries and the Chabrier IMF to model the stellar components. To model dust attenuation in galaxies ProSpect uses the Charlot & Fall (2000) model. prospect utilizes the Dale et al. (2014) templates to model the re-radiation of photons absorbed by dust into the infrared. The authors, point out that, with respect to other SED fitting code, e.g. MAGPHYS(da Cunha et al. 2008), based on the same templates, the benefit of ProSpect lies in the fact that it is extremely flexible in how it processes star formation histories and on the fact that it incorporates evolving galaxy metallicities. The author state also that such differences in the SED fitting approach lead to an increase of the stellar mass estimates by 0.2 dex. The comparison with the MAGPHYS stellar masses over the GAMA fields is shown in Section 5 of Robotham et al. (2020). No discrepancies, instead, are reported with respect to the galaxy SFR. This is also shown in the comparison provided by Thorne et al. (2021) with the MS estimated by Leslie et al. (2020) based on radio emission. These are, in turn, in agreement with the 𝑈𝑉 , 𝐼 𝑅 based MS of Whitaker et al. (2014) and Popesso et al. (2019a) over the same area. Thus, we correct only the stellar masses by 0.2 dex to bring the MS estimates of Thorne et al. (2021) on the same calibration framework of the  \nFigure A1. Comparison of the MS estimates of Leja et al. (2022) based on the SED fitting code Proppector (empty triangles) and those of Whitaker et al. (2014) based on the combination of 𝐼 𝑅 , 𝑈𝑉 in the same deep fields (filled points). The stellar masses of Leja et al. (2022) are corrected by 0.3 dex as indicated in Leja et al. (2019), for consistency with the stellar masses of Whitaker et al. (2014), based on FAST. No correction is applied to the SFR estimates. The MS are color coded as a function of the Universe age as indicated in the figure and are displaced by 0.4 dex from one another for clarity. The solid line indicates the best fit of the MS as a function of stellar mass and time based on Eq. 14, while the dashed lines indicate the best fit 1 𝜎 error.  \n<!-- image -->  \nother publications. We correct also the median to the mean and to the Kroupa IMF.  \n- GLYPH<15> Sherman et al. (2021) provide the MS of massive galaxies at stellar masses larger than 10 11 𝑀 GLYPH<12> at 1 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 0 over an area of 17.5deg 2 . They use 𝐸𝐴𝑍𝑌 GLYPH<0> 𝑃𝑌 to perform SED fitting using twelve Flexible Stellar Population Synthesis (FSPS Conroy et al. 2009; Conroy & Gunn 2010) templates in the non-negative linear combination. The 𝐸𝐴𝑍𝑌 GLYPH<0> 𝑃𝑌 FSPS templates are built with a Chabrier IMF, Kriek & Conroy (2013) dust law, solar metallicity, and star-formation histories including bursty and slowly rising models.  \nStar-forming galaxies are identified by looking at the SFR distribution per stellar mass bin. The authors identify the minimum corresponding to the green valley galaxies to isolate the star-forming systems above this minimum SFR. The MS is estimated as the average SFR of the star-forming galaxies in 4 stellar mass and 4 redshift bins. We correct to the Kroupa IMF.  \n- GLYPH<15> Leja et al. (2022) use the panchromatic SED-fitting code Prospector to measure the star-forming Main Sequence across 0 GLYPH<149> 2 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 0 using the COSMOS-2015 and 3D-HST UVIR photometric catalogs. In particular, they measure the starforming sequence using stellar population properties inferred by the Prospector𝛼 model built in the Prospector SED-fitting code. The model has 14 free parameters, consisting of a 7-component non-parametric star formation history using the 'continuity' prior which disfavors sharp changes in SFR(t) (Leja et al. 2019), a twocomponent dust attenuation model with a flexible dust attenuation  \ncurve (Noll et al. 2009), free gas-phase and stellar metallicity, and a mid-infrared AGN component with a free normalization and dust optical depth. Prospector includes dust emission powered by energy balance (da Cunha et al. 2008) with a SED and nebular emission self-consistently powered by the model stellar ionizing continuum (Byler et al. 2017). As shown in Leja et al. (2019), Prospector , or its incarnation in Prospector𝛼 , provides larger stellar masses with respect to previous SED fitting code as CIGALE (Noll et al. 2009; Boquien et al. 2019; Yang et al. 2020), MAGPHYS (da Cunha et al. 2008) or FAST (Kriek et al. 2009), largely used in the other publications collected here (see Appendix B for a discussion on the cause of the discrepancy). The authors indicate that such stellar mass increase is on average 0.3 dex with respect to other codes. Leja et al. (2022) also indicate that the SFRs derived by Prospector are lower with respect to the SFRs based on the 𝐼 𝑅 , 𝑈𝑉 contribution, derived in Whitaker et al. (2014). However, as shown in Fig. A1, we find that after correcting the stellar masses by 0.3 dex, to this quantity to the same calibration framework of all other publications, we do not observe any significant discrepancy with respect to the SFR estimated with different indicators, e.g. the 𝐼 𝑅 , 𝑈𝑉 indicator used in Whitaker et al. (2014). Indeed, the SFRs of Leja et al. (2022) are consistent with the MS best fit as a function of redshift and stellar mass, within the observed scatter of 0.08 dex. Thus, we do not apply any correction to the SFR if not to for the Kroupa IMF.  \nTheyuseaflexible neural network known as a normalizing flow to identify the ridge line of the density distribution in the 𝑙𝑜𝑔 ' 𝑆𝐹𝑅 ' GLYPH<0> 𝑙𝑜𝑔 ' 𝑀 GLYPH<3> ' plane.", 'APPENDIX B: THE KE12 SFR CALIBRATION': "KE12 offer a collection of different calibrations which leads to consistent SFRs based on different SFR indicators. In this work, we rely mostly on the calibration of Murphy et al. (2011) and Hao et al. (2011) for the 𝐼 𝑅 , 𝑈𝑉 , 𝐼 𝑅 , H 𝛼 , and radio 1.4 GHz SFR indicators. Murphy et al. (2011), in particular, uses the free-free emission measured in the Ka GLYPH<0> band (26 GLYPH<0> 40 GHz) for 10 star-forming regions in the nearby galaxy NGC 6946, including its starbursting nucleus, to compare a number of SFR diagnostics. These diagnostics include non-thermal radio (i.e. 1.4 GHz), total infrared (IR 8 GLYPH<0> 1000 𝜇 m), and warm dust (i.e., 24 𝜇 m) emission, along with hybrid indicators that attempt to account for obscured and unobscured emission from starforming regions including 𝑈𝑉 , 𝐼 𝑅 measurements. By construction, the SFR derived from the indicators calibrated in this way, result to be consistent with each other. Nevertheless, the KE12 calibration is based on the data obtained at 𝑧 GLYPH<24> 0, and one might ask if this calibration still holds for high redshift galaxies.  \nThere is substantial disagreement between the observed SFRs based on the KE12 or similar calibrations, and the SFRs estimated in simulated galaxies, as provided in state-of-the-art hydrodynamical simulations. The discrepancy is that the observed SFRs tend to be larger with respect to the predicted ones, in particular, at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3. This questions the reliability of the observed SFRs (Hayward et al. 2014; Katsianis et al. 2016; Davies et al. 2016; Martis et al. 2019; Leja et al. 2019; Lower et al. 2020). Katsianis et al. (2020), for instance, argue that the SFR based on the combination of the UV emission and the Spitzer 24 𝜇 mor Herschel 70, 100, 160, 250, 350 𝜇 m emissions are largely overestimated up to 1 dex with respect to the SFRs of simulated galaxies. This work compares the intrinsic SFRs of EAGLE simulated galaxies and those obtained through the KE12 or similar calibrations by analyzing the synthetic  \nSEDgeneratedfor the same systems by Camps et al. (2018). Furthermore, Leja et al. (2022) argue that the galaxy star formation histories reconstructed by the SED fitting code Prospector (Leja et al. 2017) lead to larger values of stellar masses and lower values of the galaxy SFRs at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3 GLYPH<149> 0, more in agreement with the SFRs predicted by simulations. This is due to the contribution of an extra component of stars older than 100 Myr, not accounted for in previous SED fitting codes. Indeed, a larger old stellar population would obviously increase a galaxy's stellar mass budget, and it might largely contribute to dust heating, thus, increasing its emission in the mid and far-infrared (see also Viaene et al. 2017; Nersesian et al. 2019; Leja et al. 2019). Leja et al. (2022) argue that neglecting such a contribution when calibrating the SFRs, derived from the 𝐼 𝑅 or 𝐼 𝑅 , 𝑈𝑉 SFR indicators, leads to a clear overestimation of the galaxy SFR as a function of redshift. For this reason, Prospector , which accounts for this contribution, provides lower SFR, bringing in agreement observed and predicted estimates at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3. However, it is worth pointing out that Leja et al. (2022) shows that such contribution is particularly relevant for galaxies below the Main Sequence, as expected (Hayward et al. 2014; Nersesian et al. 2019). It is less significant in the MS region, where the larger old stellar population affects mostly the galaxy stellar mass rather than the SFR. This is the main reason why the MS estimated by Leja et al. (2022) are in agreement with previous estimates, based on KE12 calibration, once the stellar masses are corrected for this discrepancy, as shown in Fig. A1.  \nTo further check if the KE12 calibration still holds at high redshift, we use the following approach. If there is an evolution in the fraction of IR emission due to star formation as a function of redshift, as proposed by Leja et al. (2022), the SFR derived through the 𝐼 𝑅 or 𝐼 𝑅 , 𝑈𝑉 emissions and other SFR indicators based on KE12, such as the non-thermal radio emission at 1.4 GHz and the nebular H 𝛼 emission, should no longer be consistent at high redshift. Indeed, the contribution of a larger old stellar population might affect the dust heating and the galaxy IR emission, but it can hardly affect the radio and the H 𝛼 emissions. The radio emission originating from star-forming galaxies is thought to be caused by Type II and Types Ib supernovae whose remnants are believed to accelerate most of the relativistic electrons in these galaxies. The same supernovae ionize the HII regions as well. Only stars more massive than GLYPH<24> 8 𝑀 GLYPH<12> produce Type II and Ib supernovae, and these have lifetimes of GLYPH<24> 3 GLYPH<2> 10 7 yr, while the relativistic electrons probably live GLYPH<24> 100 yr. Radio observations, therefore, probe very recent star formation in galaxies, and have the advantage that the contribution to radio emission of stellar populations older than 100 Myr is insignificant. Similarly, the H 𝛼 emission in star-forming galaxies is dominated by nebular emission in the HII regions, ionized by the most massive stars which have a short lifetime (Shivaei et al. 2015b; Katsianis et al. 2016). The strongest possible contamination is usually due to the AGNemission,ratherthanoldstars. Those might contribute through the emission generated by stellar winds, but this is estimated to be negligible as seen in the local SDSS galaxies. Indeed, it is usually seen in evolved green valley galaxies not classified as star-forming in the BPT diagram (Concas et al. 2017).  \nWe, thus, compare the SFRs derived from the 𝐼 𝑅 , 𝑈𝑉 and those derived from the radio 1.4 GHz and H 𝛼 indicators, calibrated with KE12, in high redshift galaxies. This is to check if they still provide consistent results or if the 𝐼 𝑅 , 𝑈𝑉 leads to an overestimation at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3. For this exercise, we use UV emission estimated at 2800 Å and the IR emission based on Spitzer 24 𝜇 m data provided in the 3D-HST CANDELS field (Whitaker et al. 2014). Those are  \nFigure A2. Left panel : Comparison of the 𝑈𝑉 , 𝐼 𝑅 based SFRs versus the SFRs based on H 𝛼 emission. The H 𝛼 emission is corrected for dust attenuation on the basis of an average value of the Balmer Decrement, calibrated as a function of the galaxy H 𝛼 luminosity and stellar mass in Domínguez et al. (2013). All SFR estimates are obtained according to the KE12 calibration. The points are color-coded as a function of the galaxy redshift. The red solid line shows the 1-to-1 relation and the dashed lines show the dispersion around the relation. Central panel : Same as in the left panel for a subsample of the 3D-HST galaxies where both H 𝛼 and H 𝛽 exhibit SNR ¡ 5 and the Balmer decrement is estimated individually for each system. Right panel : Comparison of the 𝑈𝑉 , 𝐼 𝑅 based SFRs versus the SFRs based on the radio emission at 1.4 GHz. The color code is the same as in the left panel. The red solid line indicates the 1-to-1 relation and the dashed lines show the dispersion around the relation.  \n<!-- image -->  \ncombined to obtain the SFRs according to Eq. 3. The stellar masses are those derived through FAST as described in Skelton et al. (2014).  \nWe use the H 𝛼 emission measured in the public 3D-HST data as described in Momcheva et al. (2016) in the same CANDELS fields. Unfortunately, the SNR of H 𝛽 , » 𝑂𝐼𝐼𝐼 … and » 𝑁𝐼𝐼 … emission lines are too low to allow any classification in the BPT diagram (Kewley et al. 2013). Thus, we select all galaxies with an H 𝛼 emission SNR higher than 5, although we are aware that there might be some contamination by AGN. The H 𝛼 emission must be corrected for dust attenuation before being converted into SFR through Eq. 4. The low SNR of the H 𝛽 emission prevents to estimate the Balmer decrement if not in a limited subsample of GLYPH<24> 80 galaxies where the SNR is higher than 5 in both H 𝛼 and H 𝛽 emission lines. Thus, we correct the H 𝛼 emission with an average estimate of the Balmer decrement as calibrated in Domínguez et al. (2013) as a function of the uncorrected H 𝛼 luminosity and stellar mass.  \nWe use the VLA-COSMOS 1.4 and 3GHz catalog (Smolčić et al. 2017) to retrieve the galaxy radio emission in the COSMOS field. This catalog provides the galaxy radio luminosity at 1.4 and 3 GHz and isolates a clean subsample of star-forming galaxies with no contamination from AGN identified in the X-rays, infrared, or through the radio excess (Delvecchio et al. 2017). The 1.4 luminosity is converted to SFR through Eq. 6.  \nThe comparison of the 𝑈𝑉 , 𝐼 𝑅 based SFRs and the SFRs derived from the radio 1.4 GHZ and H 𝛼 emissions is shown in Fig. A2. The comparison with the H 𝛼 -based SFRs is shown in the left panel. The points are color-coded as a function of the redshift, as indicated in the figure. We do not observe a clear overestimation of the KE12 𝑈𝑉 , 𝐼 𝑅 based SFRs with respect to the H 𝛼 based SFRs. All data points are distributed around the 1-to-1 relation with a scatter of 0.25 dex, without any clear bias as a function of redshift. On the contrary, we observe outliers with larger values of the H 𝛼 based SFRs with respect to the 𝑈𝑉 , 𝐼 𝑅 based SFRs, which we interpret as possible AGN contaminants. The central panel shows the same comparison for the subsample of 3D-HST galaxies where the H 𝛽 emission line SNR allows measuring the Balmer Decrement.  \nAlso in this case, we do not observe a clear overestimation of the 𝑈𝑉 , 𝐼 𝑅 -based SFRs. Nevertheless, we point out that the uncertainties are large also because we do not attempt any correction for absorption and possible line blending issues as done in Domínguez et al. (2013). The right panel shows the comparison with the SFR based on the radio emission calibrated according to KE12. All data points are distributed around the 1-to-1 relation with a scatter of 0.22 dex and without any clear over- or underestimation as a function of redshift. Few residual sources exhibit a radio excess with respect to the 𝑈𝑉 , 𝐼 𝑅 -based SFRs, which we interpret as radio AGN contaminants not included in the catalog of Delvecchio et al. (2017).  \nWe, thus, conclude that there is no clear evidence for an overestimation of the 𝑈𝑉 , 𝐼 𝑅 based SFRs at 0 GLYPH<149> 5 GLYPH<159> 𝑧 GLYPH<159> 3 and that the KE12 calibration allows obtaining consistent SFR estimates up to 𝑧 GLYPH<24> 3 from the different SFR diagnostics considered here.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author."}