The ALMA REBELS survey: obscured star formation in massive Lyman-break galaxies at z = 4–8 revealed by the IRX– 𝛽 and M ★ relations

We investigate the degree of dust obscured star formation in 49 massive (log 10 ( 𝑀 ★ / M ⊙ ) > 9) Lyman-break galaxies (LBGs) at 𝑧 = 6 . 5–8 observed as part of the ALMA Reionization Era Bright Emission Line Survey (REBELS) large program. By creating deep stacks of the photometric data and the REBELS ALMA measurements we determine the average rest-frame UV, optical and far-infrared (FIR) properties which reveal a significant fraction ( 𝑓 obs = 0 . 4–0 . 7) of obscured star formation, consistent with previous studies. From measurements of the rest-frame UV slope, we find that the brightest LBGs at these redshifts show bluer ( 𝛽 ≃ − 2 . 2) colours than expected from an extrapolation of the colour-magnitude relation found at fainter magnitudes. Assuming a modified blackbody spectral-energy distribution (SED) in the FIR (with dust temperature of 𝑇 d = 46 K and 𝛽 d = 2 . 0), we find that the REBELS sources are in agreement with the local “Calzetti-like” starburst Infrared-excess (IRX)- 𝛽 relation. By reanalysing the data available for 108 galaxies at 𝑧 ≃ 4–6 from the ALPINE ALMA large program using a consistent methodology and assumed FIR SED, we show that from 𝑧 ≃ 4–8, massive galaxies selected in the rest-frame UV have no appreciable evolution in their derived IRX– 𝛽 relation. When comparing the IRX– 𝑀 ★ relation derived from the combined ALPINE and REBELS sample to relations established at 𝑧 < 4, we find a deficit in the IRX, indicating that at 𝑧 > 4 the proportion of obscured star formation is lower by a factor of ≳ 3 at a given a 𝑀 ★ . Our IRX– 𝛽 results are in good agreement with the high-redshift predictions of simulations and semi-analytic models for 𝑧 ≃ 7 galaxies with similar stellar masses and SFRs.


INTRODUCTION
The onset of dust creation represents a milestone in the history of the Universe, as it relies on the adequate enrichment of the galaxy with metals, formation of sufficient dust particles in high-redshift supernova and inter-stellar medium (ISM) properties conducive to the survival (and growth) of dust (e.g.Draine 2003;Mancini et al. 2016;Gall & Hjorth 2018;Leśniewska & Michałowski 2019;Graziani et al. 2020;Dayal et al. 2022;Di Cesare et al. 2023).The presence of dust within galaxies can be detected through the reddening of the restframe UV and optical light in addition to emission in the mid and far-infrared (FIR).The measurement of the rest-frame FIR modified blackbody emission provides a direct signal of the presence of dust, whereas changes in the rest-frame UV colours can be attributed to other properties of the galaxy such as older ages and an increased metallicity.The majority of the highest-redshift galaxies found within deep optical to near-infrared (NIR) surveys have been shown to have ★ rebecca.bowler@manchester.ac.uk blue rest-frame UV slopes (parameterised as   ∝   ,  ≃ −2), leading to the inference of young ages and low dust content (Dunlop et al. 2012;Bouwens et al. 2014).In the past decade however, the direct detection of dust continuum emission in individual or small samples of  ≳ 7 galaxies (e.g.Tamura et al. 2019;Wong et al. 2022;Hygate et al. 2023;Hashimoto et al. 2023) has revealed the presence of dust within galaxies less than 800 Myr after the Big Bang.
Observations of star-forming galaxies in the rest-frame FIR have demonstrated the key importance of considering dust obscured starformation in galaxy evolution, with more than half of ongoing star formation being obscured at cosmic noon ( ≃ 3; see review by Madau & Dickinson 2014).There is evidence that obscured star formation continues to be important, and potentially dominates the total cosmic SFR density (CSFRD) in the range 3 <  < 6, from measurements based on rest-frame UV selected samples (e.g.Novak et al. 2017;Khusanova et al. 2021) and highly dust-obscured galaxies including serendipitous objects (e.g.Gruppioni et al. 2020;Talia et al. 2021;Loiacono et al. 2021), as well as deep ALMA and radio surveys (e.g.Zavala et al. 2021;van der Vlugt et al. 2022).Recent results extending these measurements to  ≃ 7 from Barrufet et al. (2023) have shown that dust obscured star-formation contributes at least 10 percent of the cosmic star-formation rate density, showing that it remains significant even into the Epoch of Reionization.These results have revealed a strong stellar mass dependence of the obscuration (e.g.Pannella et al. 2009Pannella et al. , 2015;;Bouwens et al. 2016;Whitaker et al. 2017), with Dunlop et al. (2017) demonstrating that at  ≃ 2 the fraction of obscured SFR rises from ≲ 0.5 at log 10 ( ★ /M ⊙ )< 9 up to 0.99 at log 10 ( ★ /M ⊙ )> 10, an effect which appears to extend to  ≃ 7 (although with a lower normalisation; Algera et al. 2023b).Direct detections of the dust continuum emission from galaxies at  > 6.5 have been made in galaxies selected from galaxies representing a wide range of intrinsic rest-frame UV luminosities for example fainter sources from lensing fields (e.g. Watson et al. 2015;Laporte et al. 2017;Tamura et al. 2019;Bakx et al. 2021;Hashimoto et al. 2023) and brighter galaxies from wide-area ground-based follow-up (e.g.Bowler et al. 2018Bowler et al. , 2022;;Schouws et al. 2022;Inami et al. 2022;Witstok et al. 2022).The obscured fraction derived from these works depends on the assumed FIR spectral-energy distribution (e.g.typically the dust temperature;  d and the emissivity index;  d ), however in general these detections reveal a significant fraction ≃ 0.2-0.8 of obscured star-formation at  ≃ 7 for galaxies of log 10 ( ★ /M ⊙ )≃ 9.5 (Dayal et al. 2022;Algera et al. 2023b).
The attenuation curve, which dictates how an intrinsic spectrum is reduced in the rest-frame UV and optical as a function of wavelength for a given optical depth, depends on the detailed properties of the dust grains and their geometric distribution (see Salim & Narayanan 2020 for a review).To directly measure the attenuation curve requires a handle on the intrinsic stellar spectra before the effect of dust, a technique that has been successfully employed at  = 2-5 (e.g.Cullen et al. 2018;Shivaei et al. 2020b).An alternative method is to compare the rest-frame UV slope, , to the ratio of the FIR to UV luminosity (infrared-excess =   = log 10 ( FIR / UV )) as the steepness of the attenuation (or extinction) curve changes the relation between IRX and  to maintain energy balance.The so called "IRX-" diagram for local starburst galaxies shows a strong correlation presented originally in Meurer et al. (1999) and then further refined in Calzetti et al. (2000).Whether this canonical Calzetti-relation holds at higher redshifts ( ≳ 2) has been the topic of debate over the past decade.An alternative to the Calzetti-like attenuation curve is the steeper relation that has been found for the Small Magellanic Cloud (SMC).Note that the SMC relation is an extinction as opposed to an attenuation curve, and there is an ongoing discussion on whether it is expected that observations of galaxies will be consistent with an SMC extinction curve when the likely complex geometry of dust is taken into account (e.g.see discussion in Cullen et al. (2018).In this case, the same column density of dust can provide an increased reddening effect in the rest-frame UV and hence a deficit from the Calzetti-like relation.Initial observations of  = 5-6 galaxies with ALMA suggested such a deficit was found (e.g.Capak et al. 2015;Barisic et al. 2017), however other studies (Bowler et al. 2018(Bowler et al. , 2022;;Schouws et al. 2022) found results consistent with the Calzetti-like IRX-.Note that the discrepancy between these studies is reduced if we consider that the first works typically assumed dust temperatures of  d = 25-45 K, while later works tended to use higher temperatures ( d ≃ 50 K).However the observation of several galaxies at  ≃ 5 that lie below the SMC prediction is still present with higher temperatures as shown by Faisst et al. 2017.Furthermore, individual sources at  ≃ 7 have been found to show significant scatter both above and below a Calzetti-like relation (Smit et al. 2018;Hashimoto et al. 2019;Bakx et al. 2020), while fainter (and likely lower mass) sources appear to show an upper limit that is even below the predic-tion of an SMC-like extinction curve (Fujimoto et al. 2016;Bouwens et al. 2016).
One key uncertainty in the measurement of the IRX- and IRX- ★ relations at increasingly high redshifts is that even at  ≳ 2 there are very few individual detections of dust continuum emission from galaxies at log 10 ( ★ /M ⊙ )< 10 (e.g.Dunlop et al. 2017;Bouwens et al. 2020).Because of this many studies at  ≳ 2 rely on stacking analyses of large numbers of individually undetected galaxies within rest-frame FIR survey data (e.g. from SCUBA-2; Koprowski et al. 2018) or alternatively small samples of often inhomogeneously detected samples from multiple follow-up programs with ALMA.In addition to the reliance on stacking or small samples, there are several systematic uncertainties that have precluded a deeper understanding of the IRX- relation at high redshift.The first is the uncertain FIR Spectral-energy distribution (SED) in the galaxies of interest, as the majority of early studies rely on measurements in the restframe FIR at typically one frequency.The derived FIR luminosity is strongly dependent on the assumed SED ( IR ∝  4+ d d ; see discussion in e.g.Behrens et al. 2018;Liang et al. 2019;Sommovigo et al. 2020) and where individual dust temperature measurements have been made a wide variation has been found in the derived  d = 20-70 K (Witstok et al. 2023).Second, the position of a galaxy with respect to the IRX- relation depends sensitively on the geometry of the dust and stars, as has been shown in theoretical works (e.g.Popping et al. 2017;Narayanan et al. 2018;Ferrara et al. 2022;Pallottini et al. 2022;Vĳayan et al. 2023).Early resolved observations have shown that there appears to be an anti-correlation between the position of the rest-frame UV and FIR emission, suggesting a complex geometry that could impact the observed IRX- (Faisst et al. 2017;Carniani et al. 2017;Bowler et al. 2018;Inami et al. 2022;Hashimoto et al. 2023;Tamura et al. 2023).
The result of these studies is an uncertain picture of how the commonly observed rest-frame UV emission in LBGs is connected to any obscured star formation at  ≳ 7 (see Hodge & da Cunha 2020 for a recent summary).To make progress in understanding obscured star formation in LBGs what is required is a statistical survey of homogeneously selected galaxies with deep observations probing the dust continuum.In this work we utilise a comprehensive survey of 49 rare, massive (log 10 ( ★ /M ⊙ )≳ 9) galaxies at  = 6.5-8.5 observed as part of the ALMA REBELS large program (Bouwens et al. 2022).The majority of these sources were selected from widearea, ground-based data over 7 deg 2 and they probe bright rest-frame UV magnitudes  UV < −21 and hence the bright-end of the restframe UV luminosity function at this epoch (e.g.Bowler et al. 2017;Harikane et al. 2022;Varadaraj et al. 2023).We also perform a consistent analysis of the ALMA ALPINE large program (Le Fèvre et al. 2020;Béthermin et al. 2020;Faisst et al. 2020b) to provide a measurement of the evolving IRX- and IRX- ★ relations from  = 4 to  = 8 using the most comprehensive homogeneous samples of  > 4 galaxies observed with ALMA.REBELS provides a unique sample of galaxies with direct dust detections (or strong upper limits) to constrain the IRX- and IRX- ★ relation within the Epoch of Reionization (EoR).This work builds upon the previous observational REBELS papers from Inami et al. (2022), Algera et al. (2023b), and Barrufet et al. (2023) and theoretical analyses tailored specifically to describe the REBELS observations from Dayal et al. (2022), Sommovigo et al. (2022a) and Ferrara et al. (2022).
The structure of this paper is as follows.In Section 2 we describe our sample from REBELS and ALPINE, presenting the ALMA observations in addition to the archival optical and NIR data that we utilize.We present the methods and results in Section 3, in particular the stacking analysis and SED fitting.In Section 4 we present the resulting colour-magnitude relation, physical properties from SED fitting and the IRX- relation.In Section 5 we discuss our results and present a new derivation of the IRX- ★ relation from  = 4-8, and we compare our ALPINE + REBELS results to the predictions from simulations in Section 6.We end with our conclusions in Section 7. Throughout this work we present magnitudes in the AB system (Oke & Gunn 1983).The standard concordance cosmology (Planck Collaboration et al. 2020) is assumed, with  0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7.

DATA
In this work we combine rest-frame UV, optical and FIR measurements to understand the dust properties of  = 4-8 galaxies.The restframe UV and optical information is provided by deep degree-scale extragalactic survey observations that have the required wavelength coverage from multiple photometric bands to select these galaxies via the redshifted Lyman break.The rest-frame FIR measurements come from targeted ALMA programs to follow-up these bright sources and provide a direct detection or upper limit on the dust-continuum emission.

REBELS
The REBELS survey is a Cycle 7 ALMA large program that observed 40 LBGs with the primary goal to measure the [CII] 158 m or [OIII] 88 m emission line.The sources were found within the ground-based Cosmological Evolution Survey (COSMOS; Scoville et al. 2007) and the XMM-Newton Large Scale Structure (XMM-LSS; Pierre et al. 2004) surveys, with the addition of two sources from HST surveys (REBELS-16 and REBELS-40).We also include 9 sources that were observed as part of the REBELS pilot programs as presented in Smit et al. (2018) and Schouws et al. (2022), resulting in a final sample containing 49 galaxies.The primary selection criterion was that the source redshift lay securely at  > 6.5 as determined by three independent SED fitting codes.The galaxies are bright in the rest-frame UV, with absolute magnitudes (measured at 1500Å in the rest-frame) in the range −23.0 <  UV < −21.3.REBELS observed each source using between two and six spectral tunings to cover the frequency range of likely FIR line emission given the photometric redshift probability distribution.As presented in Bouwens et al. (2022), Inami et al. (2022), Schouws et al. (2023, in prep), andvan Leeuwen et al. (2023, in prep), 25 galaxies have been spectroscopically confirmed via the [CII] line (with no [OIII] confirmations to-date).In addition, these observations simultaneously allowed a measurement of the dust-continuum emission, with 18 sources detected in the continuum at > 3.3 by Inami et al. (2022).The [CII] line (if detected) was masked in the continuum images.The typical continuum depth of the Band 6 or 7 data (approximately 240 GHz and 350 GHz depending on exact line scan frequencies) was  rms = 10-20 Jy with a beam of 1.2-1.6 ′′ full width at half maximum (FWHM).

ALPINE
The ALPINE survey is an ALMA large program awarded in Cycle 5 that aimed at detecting the [CII] line and dust continuum emission in 118 galaxies at  = 4-6 (Le Fèvre et al. 2020;Béthermin et al. 2020).The sources were selected from a large red-optical spectroscopic survey of "normal" star-forming galaxies within the COSMOS and Extended-Chandra Deep Field South (ECDFS) fields.The resulting sample consisted of 67 sources with spectroscopic redshifts from Lyman− emission or rest-frame UV absorption features between  = 4.4-4.6 and 51 objects at  = 5.1-5.9(Faisst et al. 2020a).Both Lyman-break selection and narrow-band selections were utilized (with more narrow-band sources in the  > 5 sub-sample), leading to a relatively high average rest-frame EW of Ly of ≃ 5-100Å (Cassata et al. 2020).In comparison to REBELS, where none of the sources had spectroscopic redshifts prior to the ALMA program, this leads to a different selection function for galaxies, which we discuss further in Section 5.As presented in Fudamoto et al. (2020b), 23 galaxies were detected at > 3.5 significance in the dust continuum in the original ALPINE survey (8 sources at  > 5).Several sources have been further followed-up in multiple bands (e.g.HZ4 and HZ6; see Faisst et al. 2020b).The typical depth of the ALPINE Band 7 data (275-373 GHz) was  rms = 30(50) Jy for the  = 5.5 (4.5) samples, with an average beam of 1.1 ′′ FWHM (Béthermin et al. 2020).

Optical and near-infrared imaging
To measure the rest-frame UV slopes of the REBELS and ALPINE galaxies, in addition to physical properties such as stellar mass, we exploited the wealth of available optical and NIR imaging data.The details of the photometry for the REBELS sample is presented in Bouwens et al. (2022) and Stefanon et al. (in prep), however we briefly describe the relevant data here.In the COSMOS (XMM-LSS) field we used the UltraVISTA (VIDEO) survey from VISTA which provided imaging in the NIR    -bands (McCracken et al. 2012;Jarvis et al. 2013).A subset of (fainter) galaxies were additionally located within a 1 deg 2 sub-region of the XMM-LSS field that has deeper observations in the  from the UK Infrared Deep Sky Survey (Lawrence et al. 2007) Ultra-Deep Survey (UDS).Spitzer/Infrared Array Camera (IRAC) photometry was extracted from the deep mosaics presented in Stefanon et al. (in prep.), in particular from the Spitzer Extended Deep Survey (SEDS; Ashby et al. 2013) and the Spitzer Matching Survey of the UltraVISTA Ultra-deep Stripes (SMUVS; Ashby et al. 2018).Photometry was extracted using 0.6 ′′ diameter apertures (0.9 ′′ for IRAC) on images where the neighbouring sources had been subtracted using MOPHONGO (Labbé et al. 2015).The aperture flux was corrected to total according to the MOPHONGO model for the galaxy.As several of the very bright ( UV < −22.5)REBELS sources were resolved in the groundbased data (Bowler et al. 2017), this step was important in deriving accurate absolute magnitudes and physical properties.Errors on the photometry were derived from empty aperture measurements on the data.
For the ALPINE sample we used the photometry provided in the COSMOS2020 catalogue (Weaver et al. 2022) which provided pointspread function (PSF) matched flux measurements for the full suite of optical and NIR filters available in the field.In comparison to the COSMOS2015 catalogue that was utilised in the original ALPINE papers (as presented in Faisst et al. 2020a), the COSMOS2020 catalogue has deeper data in a range of filters.Particularly for this work, the optical (from Subaru) and near-infrared (from UltraVISTA DR4) data are up to 1 magnitude deeper, providing significantly improved measurements of the rest-frame UV slope and derived stellar masses.We used the 'Classic' catalogue that provides photometry measured in 2 ′′ diameter apertures.In the COSMOS2020 catalogue this aperture photometry is corrected to a total flux by a constant factor determined from the Source Extractor MAG_AUTO.From the full ALPINE catalogue of 118 sources, ten sources lie within The average optical and NIR photometry and corresponding bestfitting SED model for the REBELS samples at  < 6.9 (top) and  > 6.9 (bottom).In each plot we show the observed fluxes and errors from the bright and faint stacks as the blue and red filled circles respectively.The six points correspond to the VISTA     and Spitzer/IRAC [3.6𝜇m] and [4.5𝜇m] bands.The open circles correspond to the synthetic photometry as derived from the best-fitting SED model from BAGPIPES, shown as the solid line.The observed change in the [3.6m] to [4.5m] colour with redshift is due to the transit of the strong H+ [OIII] emission lines between these filters (e.g.Smit et al. 2014).
the ECDFS.To provide a sample with uniform photometry from the COSMOS2020 catalogue, we excluded these ten sources from further analysis, leaving a final ALPINE sample of 108 galaxies.

METHODS
The primary goal of this work is to measure the rest-frame UV and FIR properties of the 49 bright  = 6-8 LBGs observed as part of the REBELS survey.We also include for comparison a consistent analysis of the ALPINE sample in the COSMOS field, to provide a base-line reaching down to  ≃ 4. We measure the properties of individual sources from both surveys, but also perform a stacking analysis to derive average properties within bins of  ★ and  UV .Due to the relatively low fraction of sources that are directly detected in the dust continuum (0.43 for REBELS, 0.19 for ALPINE), stacking is a key tool to understand the dust continuum properties of galaxies within this sample.Here we describe the key methods used in this work.Where possible we used identical approaches for the REBELS and ALPINE analysis to provide a direct comparison between galaxies spanning  ≃ 4 to  ≃ 8.

Stacking analysis
The ALPINE and REBELS samples span a range of redshifts leading to different rest-frame features being observed in our available optical and NIR data.We therefore separate our sample into two redshift bins in ALPINE and three within REBELS.The two ALPINE bins correspond to  = 4.0-4.5 and  = 4.8-5.4,which are the two main redshift groupings within the sample.All ALPINE sources were observed in the same ALMA band (Band 7).In REBELS we split the sample into three bins, by increasing redshift as detailed below.We excluded the HST selected sources REBELS-16 and REBELS-40 from the stacks as they have different rest-frame optical and NIR filters available.The first bin included galaxies at  = 6.5-6.9 (20 galaxies), with the second bin included galaxies in the range  = 6.9-7.7 (20 galaxies).This separation at  = 6.9 was chosen as it is the point at which the Lyman break starts to move into the VISTA  -band and is also when the H+[OIII] 4959, 5007 rest-frame optical lines move from the [3.6m] to the [4.5m] band.The 40 sources within these two bins had observations in ALMA Band 6.The third and final REBELS bin included the seven galaxies that have  > 7.7.Four of these sources have ALMA Band 7 observations that were designed to target the [OIII] line.These sources have photometric redshifts in the range  = 7.7-8.6,however none have been spectroscopically confirmed to-date.Due to the wide range of photometric redshifts in this bin, and the differing ALMA measurement bands, we do not create a stack from this sub-sample.However, we present their individual IRX- properties for comparison with our stacks.We further choose to split the sample using  UV as this provided the greatest dynamic range in the derivation of the IRX- and IRX- ★ relation, while also not suffering from biases (as the  and  ★ values have significant statistical errors, scatter between bins can lead to biases; see e.g.McLure et al. 2018).The bins we used for REBELS are shown in Table 1.We split the sub-samples by  UV = −22.0( UV = −22.5)at  ≃ 6.7 ( ≃ 7.2) to provide roughly equal sources with the brighter/fainter absolute magnitude bins.For ALPINE we split by  UV = −22.0for both redshift bins, and in addition we separate into two stellar mass bins to take into account the wider range of  ★ values within the sample.The ALPINE bins are detailed in Table A1.By restricting the  UV and  ★ ranges slightly, we result in a final ALPINE sample of 54 galaxies within the  ≃ 4.5 bin and 32 within the  ≃ 5.5 bin.
We performed weighted mean stacking of the ALMA continuum data in the specified bins in the image plane.Our results are unchanged if we instead use a median stacking procedure.Due to the majority of sources in the REBELS and ALPINE samples being undetected in the ALMA data we stack at the position of the observed rest-frame UV.We note that this could cause an underestimate of the peak ALMA flux if there exists significant offsets between the restframe UV and FIR emission (e.g. as simulated in Bowler et al. 2018).Given the large beam of the REBELS and ALPINE observations, and the relatively small absolute offsets found for these samples (e.g.Le Fèvre et al. 2020;Inami et al. 2022, we expect the effect to be small.Indeed, as described in Section 3.4, only the most massive ALPINE stack shows evidence for extended emission, which we attribute partially to offsets between the rest-frame UV and FIR positions.We take this into account in the flux measurement by using a Gaussian fit to the stack. Foreground sources were masked based on a map of objects derived from PyBDSF (Mohan & Rafferty 2015), excluding pixels within a radius of 1.5 ′′ from the rest-frame UV centroid of the galaxy.The majority of sources did not have any neighbours at the depth of the ALMA imaging.Errors on the stacked flux measurements were determined by remaking each stack using bootstrap resampling with replacement.To determine the mean optical and NIR photometric measurements we averaged the fluxes from the REBELS and COSMOS2020 catalogues (Section 2.3).We also experimented with stacking the images themselves, however we concluded that it was not possible to improve on the flux stacking results, due to close neighbours contaminating the flux measurements.This contamination was taken into account in our catalogue creation, with nearby sources being subtracted prior to aperture photometry (see Section 2.3).For ALPINE we use the COSMOS2020 'Classic' aperture photometry measurements where basic subtraction of neighbours is performed.We visually checked the ALPINE sources in the COSMOS optical and NIR imaging, but found that they are all sufficiently isolated for the catalogue fluxes to be robust.

SED fitting
We fit the photometric data for the individual REBELS and ALPINE sources (and the derived stacks) using BAGPIPES (Carnall et al. 2018) to provide a best-fitting model with which to measure the restframe UV slope.The fitting also provides physical properties for the stacks, which we include in particular for measuring the IRX- ★ relation.We fix the redshift to the spectroscopic redshift when available (28 sources in REBELS, all of the sources in ALPINE), and for the stacked photometry we fix the redshift to the average redshift.We found that using a luminosity-weighted redshift instead of an average had no effect on our results, as the difference was  ≤ 0.03.We include bands above the Lyman-break reaching to the [4.5m] filter, beyond which the resolution and depth decreases dramatically.We also exclude bands that contain the Lyman-break in the fitting of the stacked photometry, as the small differences in break position within the band lead to tensions within the fitting.Hence for REBELS, we fit to the    [3.6m] [4.5m] bands for the  = 6.5-6.9 sub-sample stack, and to the   [3.6m] [4.5m] bands for the  > 6.9 stack.The resulting photometry and best-fitting SED models for the REBELS stacks are shown in  2022) can increase the derived stellar masses by on average ≃ 0.5 dex and in some cases ≳ 1.5 dex.To provide a closer comparison to previous literature measurements of the IRX- ★ we primarily consider the  ★ values derived with standard parametric 6.5 < z < 6.9 23.5 < M UV < 22 6.5 < z < 6.9 22 < M UV < 20.5 6.9 < z < 7.7 23.5 < M UV < 22.5 6.9 < z < 7.7 22.5 < M UV < 20.5 SFHs, however we note where relevant how our results would change with an assumed alternate SFH.For the ALPINE sample, Faisst et al. (2020b) found that the -value derived depended on the assumed dust law in the fitting.We also recover this trend in our sample, with the dervied -slopes being redder by around 0.1 when fitting with an SMC dust law in comparison to a Calzetti law.

Rest-frame UV luminosity and slope determination
The monochromatic rest-frame UV luminosity was derived at 1500 Å using a top-hat filter of width 100 Å applied to the best-fitting SED model from BAGPIPES for both the REBELS and ALPINE stacked photometry.We note that the aperture photometry for both the REBELS catalogue and COSMOS2020 have been corrected to a total flux accounting for the full extent of the galaxy, and hence the derived  UV can be considered a total absolute magnitude.We found a systematic offset brightwards of Δ UV = −0.1 mag between the ALPINE absolute magnitudes presented in Faisst et al. (2020b) and those found in our analysis, with some sources having considerable offsets reaching > 0.5mag.Further inspection reveals this to be due to the improved photometry between the COSMOS2015 and COS-MOS2020 catalogues.The rest-frame UV slope is historically defined from a series of windows in the continuum from  rest = 1268-2580 Å (Calzetti et al. 1994).Different methods for measuring  from the available photometric data in high-redshift galaxies have been extensively discussed, including the fitting of a power law or a Table 1.The measured FIR fluxes and derived properties of the four REBELS stacks at  = 6.5-7.7.The equivalent results for the ALPINE analysis is presented in Table A1.The top (bottom) two rows show the lower (higher) redshift stack, with the stacks ordered by  UV .Columns 1 and 2 detail the redshift and number of sources included in each stack.The average redshift and  UV of each stack are shown in Columns 3 and 4. Column 5 presents the measured peak ALMA flux, with the corresponding S/N shown in brackets.The flux measured using a Gaussian fit is shown in Column 6.The derived FIR luminosity (assuming  d = 46 K,  d = 2.0) and the resulting IRX value are shown in Columns 7 and 8.The  IR was determined from the peak flux for the REBELS results.Finally the rest-frame UV slope  is presented in Column 9, as measured from the best-fitting SED model.
/ Jy / Jy /10 11 L ⊙ 6.5 < z < 6.9 7 6.76 −22.34 ± 0.14 30.3 ± 7.9 (5.6) 31.1 ± 9.9 1.4 ± 0.4 −0.11 +0.16   −0.18 −2.10 +0.14 −0.14 6.5 < z < 6.9 13 6.67 −21.71 ± 0.12 44.5 ± 10.4 (6.5) 64.0 ± 15.3 2.1 ± 0.5 0.31 power law with a Lyman-break to the photometry directly or to the fit SED model (e.g.Dunlop et al. 2012;Rogers et al. 2013).In this work we measure  from the best-fitting SED model derived from BAGPIPES (see Section 3.2), excluding regions that are outside of the Calzetti windows, to avoid strong absorption or emission features.Errors were derived using a bootstrap analysis, where we restacked the photometry and re-fit using BAGPIPES.Comparing the derived  values to those presented in the original ALPINE analysis (Faisst et al. 2020b), with find on average a very mild bias to bluer slopes by 0.05 in our analysis, when comparing results obtained by fitting assuming the same dust attenuation law.The scatter between individual objects can be large (up to  = 0.5), but is within the errors of the derived  values.

Rest-frame FIR luminosity derivation
Using PyBDSF (Mohan & Rafferty 2015) we measured both the peak flux and flux derived from a Gaussian fit for the individual sources and the stacked ALMA data.We found that our stacked results were consistent with being unresolved for log 10 ( ★ /M ⊙ )< 10 (i.e. the full REBELS sample and low-mass sub-sample of ALPINE).We define an unresolved source if the measured major and minor axes from PyBDSF are consistent with the beam size within the 1 error (again derived within PyBDSF).Hence we used the peak flux measure in these cases.For the ALPINE sample at log 10 ( ★ /M ⊙ )> 10 we instead used the Gaussian flux measurement, as the derived sizes from PyBDSF were significantly resolved.From the single data-point in the observed mm-regime from ALMA in Band 6 for ALPINE and Band 6 or 7 for REBELS we determined the total FIR luminosity by assuming a modified blackbody SED.We corrected for the effect of the Cosmic Microwave Background following da Cunha et al. ( 2013), which results in an increase in the  IR by 10 percent for the REBELS sample.While some sources within the two surveys have been observed in multiple ALMA bands (Algera et al. 2023a), we choose here to provide a uniform measure of  IR from the single main band that is available for the full REBELS + ALPINE samples.
For the analysis presented in this work we assumed a single fixed dust temperature of  = 46 K, with an opacity fixed to  d = 2.0 (consistent with that recently measured by Witstok et al. 2023).This dust temperature was derived by the model of Sommovigo et al. (2022a) and was used by Inami et al. (2022)

RESULTS
In Fig. 1 we present the stacked photometry and best-fitting SED model for the REBELS sample, split into the two main redshift (and further two  UV bins) as shown in Table 1.The results of stacking the ALMA data in these bins are shown in Fig. 2. We find a significant (7-10) detection in the fainter  UV -bin for both the  ≃ 6.7 and the  ≃ 7.2 stacks.In the brighter stacks we find marginal detections in the dust continuum, at 4 and 2.5 for the  ≃ 6.7 and the  ≃ 7.2 stack respectively.The fluxes we derive are consistent with that found in the independent analysis of the REBELS sample by Algera et al. (2023b), who used a Monte Carlo stacking analysis to measure a correlation of  IR with stellar mass.From these ALMA detections we then proceeded to compute the  IR and combine this with the restframe UV information ( UV , rest-frame UV slope) and the stellar mass as derived from BAGPIPES as detailed below.

Physical properties
The splitting of the bright REBELS sample into two redshift bins separated at  = 6.9 allows us to provide high-S/N stacks of the rest-frame UV and optical emission where the strong

Colour-magnitude relation at z ≃ 7
In Fig. 3 we present the rest-frame UV slopes of the REBELS sources, and the stacks, in comparison to colour-magnitude relations from recent literature measurements.The REBELS sample allows us to measure the colour-magnitude relation up to  UV ≃ −23, which is considerably brighter than the typical galaxies found and studied previously using HST or JWST data that are typically dominated by sources at  UV ≳ −21.We find that the REBELS galaxies show a range of rest-frame UV slopes, with −2.7 <  < −1.0, with many measured -slopes dominated by large errors (Δ > 0.5).Reassuringly, our  values derived from the stacked photometry follow the distribution of individual values.We compare our results to the extrapolated relations from fainter sources derived at  ≃ 7 using HST data by Bouwens et al. (2014), and the two recent JWST results by Topping et al. (2023) at  = 7.3 and by Cullen et al. (2023) at  = 8-12 (here we use the redshift evolution applied to the slope found Rogers et al. 2014).These studies computed  by fitting a power law the available photometry probing the rest-frame UV.At  ≃ 6.7 we find good agreement with these relations in our fainter bin, however we see that our brighter stack has a significantly bluer slope than expected from the extrapolated colour-magnitude relations from previous studies at fainter magnitudes.At  UV < −22 the offset bluewards from the colour-magnitude relations is between Δ ≃ 0.3-0.7 depending on the study.Looking at the slightly higher redshift bin at  ≃ 7.2 we find an offset to bluer  values by Δ = 0.4 in both the brighter and fainter stack.Our results support a flattening, and potentially even a turn-over, of the colour-magnitude relationship at  UV ≲ −22, with these galaxies showing a mean colour of  = −2.1 in contrast to the predicted colour of  ≃ −1.4 to −1.7 from the relations extrapolated from fainter LBGs.As we discuss further in Section 5, this turn over can be explained by the effect of scatter in the obscuration when considering sources that have a steeply declining number density.
In the measurement of the colour-magnitude relationship we must consider any effect of sample selection and -measurement bias in the results we obtain.It is possible that we could be missing redder  ≃ 7 galaxies due to the requirement that the sources show good high-redshift fits and poorer quality fits to (typically redder) lowredshift galaxy contaminants.As shown in Fig. 3 we are able to measure -values as red as  ≃ −1.2 for the sources in our sample, even at the faint end, whereas we do not find significant numbers The REBELS galaxies at  = 6.5-6.9 (upper) and 6.9 <  < 7.7 (lower) in comparison to the colour-magnitude relation found by previous studies.In each plot the individual galaxy measurements of the rest-frame UV slope and  UV are shown as the grey open points, while the stacked results (from the fits shown in Fig. 1) are shown as the blue filled points.The relationship derived from fainter studies are shown as the black solid, dotted and dashed lines from Cullen et al. (2023), Topping et al. (2023) and Bouwens et al. (2014) respectively.Slight differences between the relations in the upper and lower plot are due to the derived evolution of the relation from these studies.We also show the data points from the same three works (at a fixed rather than evolving redshift), note in particular that the redshift range of the data from the Cullen et al. ( 2023) study is at  > 8.
of the brightest sources to be as red (even though the increased S/N should make bright, red, sources easier to identify than similarly red, fainter sources).The REBELS sample selection is not only based on the rest-frame UV bands but also includes the [3.6m] and [4.5m] bands in the SED fitting.The Spitzer/IRAC colour aids in the selection of robust  ≃ 7 galaxies due to the specific colours produced by the rest-frame optical nebular emission lines in the [3.6m] and [4.5m] bands (Smit et al. 2015).Using these filters could be biasing our sample towards bluer slopes by potentially selecting young galaxies with stronger nebular emission.We discount this however, as the distribution of the EW 0 (H + [OIII]) of the REBELS sample is in excellent agreement with  ≃ 7 samples that are selected only based on a strong Lyman break (see figure 18 Table 2.The physical properties of the REBELS sources as derived from the stacked photometry.We employ both SFR calibrations and SED fitting using BAGPIPES, where the best-fitting models are shown in Fig. 1.Each row corresponds to a different stack in redshift and absolute UV magnitude as shown in Columns 1 and 2 respectively.In Columns 3 and 4 we present the SFR derived from the rest-frame UV and FIR (see Section 4.1), with Column 5 showing the obscured SFR fraction derived from these quantities.Columns 6 to 10 show the SFR,  ★ , age,  V and ionization parameter as derived from BAGPIPES assuming a SFH following a delayed  model and a fixed metallicity of  = 0.2 Z ⊙ .Non-detections Individual detections REBELS (6.5 < z < 6.9) REBELS (6.9 < z < 7.7) of Bouwens et al. 2022).In fact, because these colours are challenging to reproduce by low-redshift galaxy contaminants it can aid in the recovery of good high-redshift galaxy fits to sources with redder restframe UV slopes (e.g. in Endsley et al. 2021; see Stefanon in prep.for individual SED fits).Hence we conclude that our measurements of the rest-frame UV slope of the REBELS LBGs are unlikely to be significantly biased, with the caveat that we only select sources that are bright in the rest-frame UV, and hence will be incomplete to the most obscured galaxies (with the extreme situation being fully 'UV-dark' galaxies as found in e.g.Fudamoto et al. 2021).

IRX-𝛽 relation at z ≃ 7 from REBELS
In Fig. 4 we present the IRX- relation derived from the REBELS sample in the redshift range 6.5 <  < 7.7.We also show the derived values for individual sources, of which there are 18 detections at > 3.3 from Inami et al. (2022).These results (both individual galaxies and for the stacks) were computed with the assumption of a modified blackbody FIR SED, with an assumed dust temperature of  = 46 K and  d = 2.0 (Section 3.4).The individual results show a large scatter horizontally on the plot as a result of the large errors in individual measurements of the rest-frame UV slope.We find that the majority of this range in observed rest-frame UV slopes in the sample can be explained with statistical scatter, with the intrinsic variation as a function of  UV derived to be of the order of Δ = 0.1-0.3 for REBELS (see Table 1).Rather the scatter can be explained simply due to the large errors on the individual  measurements, which we have demonstrated via a simple simulation assuming the sample is drawn from a constant input  = −2.0 with the same  measurement errors.With our assumed rest-frame FIR SED, based on the work of Sommovigo et al. (2022a), we do not confirm any sources significantly below the SMC-like IRX- relation as found by previous high-redshift studies (e.g.Barisic et al. 2017;Faisst et al. 2017;Smit et al. 2018; note that these works assumed a lower  d ≃ 30-40).Although at the depths of our observations, 31 of the 49 sources in REBELS are undetected in the dust continuum and hence the IRX values represent upper limits.We see one source that is significantly in excess of the others, with IRX = 1.2 ± 0.2.This is the unusually FIR bright object REBELS-25 that is discussed further in Hygate et al. (2023).This source is included in our stacks, however our results are unchanged if it is removed.These data are shown in comparison to the expected relation for a Calzetti et al. (2000) dust attenuation and SMC dust extinction law, assuming an intrinsic -slope,  0 = −2.3.This intrinsic slope is consistent with that found in detailed SED fitting of comparable mass sources at  = 3 (McLure et al. 2018) and similar to that found in simulations of galaxies at  ≃ 5 (e.g.Cullen et al. 2018 found  0 = −2.4).We present a fit to the IRX- relation, and fits from previous works at higher redshift, that include a steeper intrinsic  in Section 4.4.We find no strong correlation between the offset from a Calzetti-like relation and  V ,  ★ , or spatial offset between the rest-frame UV and FIR flux (from Inami et al. 2022).There is a weak trend that the FIR brightest galaxies tend to be above the relation, such as REBELS-25 (as has been seen for ULIRGS; see discussion below).
Turning to the stacked results, we present four individual points corresponding to the two different redshift and  UV bins.For the 6.5 <  < 6.9 sub-sample we see that the brighter and fainter stack shows significantly different rest-frame UV slopes, with the C a lz e t t i S M C REBELS (6.5 < z < 6.9) REBELS (6.9 < z < 7.7) Bowler+22 COS-87259 (Endsley+23) Witstok+22 A1703-zD1 (Molyneux+22) A1689-zD1 (Bakx+21) Schouws+22 brighter stack ( UV < −22) appearing bluer, while simultaneously being fainter in the FIR.The brighter stack also shows a lower ALMA detected flux, and hence a lower IRX both from a higher  UV and a reduced  IR .The same trend is seen for the galaxies in the 6.9 <  < 7.7 sub-sample, however here as both stacks are blue in the rest-frame UV (and the brighter stack is undetected in the FIR), we have a reduced dynamic range in .Overall, we find, somewhat counter-intuitively to the consensus colour-magnitude relation (Section 4.2 and Fig. 3), that the rest-frame UV brightest galaxies in REBELS are bluer than the sources at slightly fainter magnitudes.As expected by the canonical IRX- relation, we find the bluer sources show a reduced IRX, and this is driven primarily due to a reduced  IR (although as we have previously described, note that galaxy age, dust SED and star-dust geometry can alter the expected relation; e.g.Popping et al. 2017).

Comparison to previous studies at 𝑧 ≃ 7
In Fig. 5 we compare our REBELS results with those from previous studies at  ≃ 7. We find that our stacked results are in good agreement with the previous measurements derived from luminous LBGs in Schouws et al. (2022) and Bowler et al. (2022).There are a handful of sources with redder rest-frame UV slopes and low IRX-values found within the study of Witstok et al. (2022), although we note that the measured -slopes are relatively uncertain in these cases.Molyneux et al. (2022) found a red rest-frame UV slope and an upper limit on the dust continuum emission in the  = 6.8 galaxy A1703-zD1.These studies all assumed a dust temperature of 50 K and  d = 1.5-1.6 and hence we expect no appreciable offset to the results of this work due to differences in the chosen rest-frame FIR SED.As shown in Fig. 8 and further discussed in Section 4.4, although the dust temperature in these cases is higher than we assume in this work, the lower  d compensates almost exactly.We additionally show two galaxies where there is a confirmed ALMA detection and robust rest-frame UV slope determination.The recent study of the  = 7.13 lensed galaxy A1689-zD1 from Bakx et al. (2021) found an IRX value that is in excess of he majority of the other points and the canonical Calzetti-like relation.The  IR of this source was derived with the observed best-fitting dust temperature of  d = 40 K, with a fixed  d = 2.03.If a higher dust temperature was assumed (to make the FIR analysis consistent with this work, and the other studies shown in this plot) this would increase the IRX by 0.2 dex (Fig. 8) resulting in an even greater excess.To show this source on the IRX- relation we take the rest-frame UV colour derived by Watson et al. (2015).Knudsen et al. (2017) argue that A1689-zD1 could be a massive starburst due to the observed [CII] deficit, large  IR (given the stellar mass) and disturbed morphology.The fact that A1689-zD1, as well as REBELS-25 in Fig. 4, appear in the upper left region of the IRX- diagram could be due to a spatial offset between the regions emitting in the rest-frame UV and FIR.The FIR emission in this case would be dominated by optically thick emission, while the rest-frame UV colour is measured from unobscured stars leading to an unusually blue colour for the observed IRX (as seen in ULIRGS; Casey et al. (2014), and predicted in theoretical works e.g.Popping et al. 2017;Behrens et al. 2018;Liang et al. 2019;Sommovigo et al. 2020;Ferrara et al. 2022).We also show the galaxy COS-87259 from Endsley et al. (2023) that was found within the COSMOS field using an LBG selection, but has been confirmed to be a highly star forming and dust obscured radio-loud AGN at  = 6.853.This source is very red, but it has a high derived IRX placing it slightly above the prediction of a Calzetti-like IRX- relation.

Individual REBELS galaxies at 𝑧 > 7.7
In Fig. 6 we show the IRX- results we derive for the seven galaxies in REBELS that have photometric redshifts at  > 7.7.These galaxies were not included in our stacking analysis due to four sources having Band 7 observations, and the relatively uncertain photometric redshifts derived for these sources at the very high-redshift end of the REBELS sample.We compare to the Hashimoto et al. (2023) work that spectroscopically confirmed a group of galaxies (nicknamed RIOJA) at  = 7.88, with three of the components showing detections in the rest-frame FIR from ALMA.Other  ≳ 7.5 sources have been observed with ALMA (e.g.MACS0416_Y1 at  = 8.31 and MACS0416-JD at  = 9.11; Hashimoto et al. 2018;Bakx et al. 2020) however these other objects do not have published rest-frame UV slopes.Two of the REBELS sources at  > 7.7 are detected in the dust continuum (REBELS-4 and REBELS-37; also called XMM-355 and UVISTA-1212 respectively in Bowler et al. 2020), while the other five are not.REBELS-4 is in good agreement with our stacked results at slightly lower redshift, however REBELS-37 shows a redder rest-frame UV colour and deficit in IRX from both the Calzetti and SMC relations shown.We note that the -slope value is more uncertain at these redshifts, due to the few bands (,   ) available for fitting, and the broader uncertainty in photometric redshift leading to a degeneracy between redshift and slope.Excluding REBELS-37, we find good agreement within the (large) errors with our  ≃ 7 stacks, although we note that the majority are upper limits on the dust continuum.C a lz e t t i S M C REBELS (z > 7.7) REBELS (6.5 < z < 6.9) REBELS (6.9 < z < 7.7) YD1,YD4,YD7 (Hashimoto+23) Figure 6.The individual IRX- points for the seven REBELS galaxies with photometric redshifts at  > 7.7 (open grey diamonds) in comparison to our stacked results at  = 6.5-6.9 (navy circles) and  = 6.9-7.7 (blue squares).We also compare to the  = 7.88 galaxy group found within the Abell cluster A2744 lensing field, that have been detected in the dust continuum by Hashimoto et al. (2023).The Calzetti and SMC-like relations are shown as described in the caption of Fig. 4.

IRX-𝛽 from z = 4 − 8 from ALPINE and REBELS
To provide a consistent comparison to the  > 6.5 results from REBELS, we performed the same analysis on the ALPINE sample in the COSMOS field (see Section 2).We present the ALPINE results compared to the individual object measurements in Fig. A1 with the values presented in Table A1.As the ALPINE sample is both larger and has a broader range of measured stellar masses, we additionally binned in  ★ .In our computation of the  IR and hence IRX for the two samples we assumed the same rest-frame FIR SED, with a dust temperature of  d = 46 K and emissivity of  d = 2.0 following the work of Sommovigo et al. (2022a,b) (see Section 3.4).We compare the ALPINE results with those from REBELS in Fig. 7.As was found for the REBELS sources, the brighter rest-frame UV stacks have lower IRX and appear bluer.We also see a strong dependence on stellar mass, with the galaxies at log 10 ( ★ /M ⊙ )> 10 showing considerably redder colours with  ≳ −1.75 and a higher IRX by ≃ 0.75.We note here that due to the selection methodology of the initial ALPINE sample, a larger fraction of the sources in the  ≃ 5.5 bin were selected to be Lyman- emitters.Taking a rest-frame EW of > 50(25) Å as the separation between LBGs and LAEs, 8 (30) percent of the  = 4.5 sub-sample are LAEs in comparison to 38 (57) percent of the  = 5.5 sub-sample.As LAEs have in general been found to show lower dust attenuation (e.g.Schaerer et al. 2015) this could explain the small offset we see between these two redshift bins.
As we discuss further in the next section, we find consistent results between the derived IRX- relation between the REBELS and ALPINE samples in our analysis, when we use the same modified blackbody SED fitting analysis and  measurement procedure.The data appears to agree with the local starburst relation of Calzetti et al. (2000), with no evidence from our stacked results for a deficit in the relation that could be consistent with SMC-like dust, given our assumptions on FIR SED.
By combining the results from the two surveys we are able to measure the IRX- relationship across a wide redshift range from  = 4-8 from the largest sample of log 10 ( ★ /M ⊙ )> 9 galaxies available with deep ALMA follow-up.Taking the stacked detections for ALPINE and REBELS, we fit the slope of the IRX- relation with a given intrinsic rest-frame UV slope,  0 , according to the formalism presented in McLure et al. (2018) as: In this formalism, the Calzetti (SMC)-like relation has a slope of d 1600 /d = 1.97(0.91)and the 1.71 pre-factor is a constant set by the bolometric correction between the total rest-frame UV emission available to heat the dust and that characterised by  UV .The prefactor can change if we break the assumption of a dust screen, however for this analysis we keep it constant.For our combined ALPINE and REBELS results we find a best-fitting slope of d 1600 /d = 2.11 ± 0.13 when assuming  0 = −2.3 or a shallower slope of the IRX- relation of d 1600 /d = 1.38 ± 0.09 when assuming  0 = −2.5.The intrinsic rest-frame UV slope of our sample is not known, however from BAGPIPES SED fitting analysis we find it to be between  0 = −2.3 and −2.5 and hence present the results of both fits.As can be seen in Fig. 7, both  0 assumptions provide a good description of the data over a broad range in measured rest-frame UV slope.Our results are in general in excess of the previously derived IRX- relations at  > 4 (e.g. from Fudamoto et al. 2020b;Schouws et al. 2022).

Comparison to previous results from the ALPINE survey
Our conclusions on the IRX- relation at  = 4-8 are different to those found in the previous ALPINE analysis presented in Fudamoto et al. (2020b), particularly at log 10 ( ★ /M ⊙ )> 10 where we find a higher IRX by around 0.5 dex when comparing stacks across the same  ★ range.The later studies of Burgarella et al. (2022) and Boquien et al. (2022) found similar conclusions to the Fudamoto et al. (2020b) study with further analysis of subsets of the ALPINE sample.Fudamoto et al. (2020b) present stacked IRX- relations using bins in  and  ★ , finding the results to be consistent.In the following discussion we compare to the  ★ binning results of Fudamoto et al. (2020b) as this has been shown to be the least biased estimator of IRX- (e.g.McLure et al. 2018).This provides the most natural comparison as our points are already stacked in  ★ , however we additionally stack in  UV bins.Hence in the following discussion we combined our  UV bins at a given  ★ .This leads to points that lie mid-way between the two  UV bins per  ★ bin, as expected.To identify the cause of this offset we first directly compared the derived -slopes,  UV values and ALMA fluxes for individual objects.We find that our rest-frame UV slopes are on average bluer than those derived in ALPINE by Δ = 0.1, with around half of this difference attributed to the dust law that we assume in the SED fitting (we assume Calzetti, whereas in Faisst et al. 2020b took the average between results with an SMC and Calzetti dust law).The  UV values are found to be offset slightly brighter (0.1 mag) in our analysis, which used the COSMOS2020 catalogue instead of the COSMOS2015 data analysed in Faisst et al. (2020b), however this has a negligible effect on the derived IRX.For the 20 percent of the ALPINE sample that have dust continuum detections we find good agreement between our raw flux measurements.Both our study and that of Fudamoto et al. (2020b) take into account the fact that the dust continuum emission may be extended in the higher mass (log 10 ( ★ /M ⊙ )> 10) stacks by using a Gaussian fit to the ALMA data.On closer inspection, the extension found in these stacked images is due to both an intrinsic extension (i.e. higher mass sources have an extended dust distribution) and an artificial extension introduced in the stacking process due to offsets between the rest-frame UV and FIR centroid.In the binning analysis, we determine the  slopes from SED fitting to the stacked optical/NIR photometry, while Fudamoto et al. 2020b take the median  in each bin.However despite the different method, when comparing the same bins in  ★ we find only a 0.1 difference between the resulting -slopes (0.2 for the lower mass bin at  ≃ 4.5), and 0.05 of the difference can be accounted for by the different assumed dust law in the fitting (Section 3).
Assuming that the fluxes in the stacks are consistent, the main difference is in the FIR SED assumed in the derivation of the  IR .Fudamoto et al. (2020b) used a scaling factor to compute  IR that was derived from an empirical FIR template created by stacking Herschel data in the COSMOS field.An expanded sample of photometrically selected galaxies over a similar redshift to the ALPINE sample was used in the creation of this template (Béthermin et al. 2020), and it can be approximated with a modified blackbody with a fixed  d = 1.8 of temperature  d = 41 ± 1 K and  d = 43 ± 5 K at  = 4-5 and  = 5-6 respectively.While only 3-5 K lower than that assumed in this work, the difference is enough to account for a 0.25 dex difference in the resulting IRX given the same input flux measurement.We mark this offset as an arrow in Fig. 7.This is illustrated in Fig. 8, where we show the offset in IRX expected for changes in  d and  d .This difference, and the slightly bluer rest-frame UV slopes we find, can account for 0.3 dex of the observed difference between our analysis of ALPINE and that presented previously in Fudamoto et al. (2020b).As can be seen in Fig. A1, the ALPINE sample consisted of 80 percent upper limits on the dust continuum emission and hence the exact stacking process could contribute to the remaining difference.Despite not knowing exactly the dust SED for the ALPINE and REBELS samples, we have shown that the two samples have consistent IRX relations when the ALMA measurements are fit with the same modified blackbody assumptions.If the ALPINE sample does indeed show a lower dust temperature to that of REBELS (as expected from the  d -redshift relation in e.g.Schreiber et al. 2018), then we would recover a lower IRX- relation for the ALPINE dataset by 0.25 dex.

DISCUSSION
This work presents a consistent analysis of the two most substantial ALMA surveys measuring the dust continuum emission from  > 4 LBGs.We have computed the IRX- relation at  ≃ 7 from REBELS, and compare this to a consistent analysis of the  = 4.5 and  = 5.5 samples from ALPINE.These surveys targeted galaxies with stellar masses of log 10 ( ★ /M ⊙ )> 9 and hence probe the high-mass end of the known galaxy distribution at these redshifts.Even before considering the results from the ALMA observations themselves, we find that the REBELS galaxies are surprisingly blue in the restframe UV as compared to the predicted colour of galaxies from the extrapolated colour-magnitude relation found for fainter sources (Fig. 3).Despite being some of the most massive galaxies known at  ≃ 7, the REBELS sources show a flatting, or a potential turnover, in the colour-magnitude relation for galaxies at  UV < −21.5.As discussed in Section 4.2, we do not think this behaviour is due to our sample selection in part due to the fact that predominantly blue slopes are found even in the most luminous (and hence highest S/N) galaxies where redder slopes would be robustly measured.Following the arguments presented in e.g.Bowler et al. 2015Bowler et al. , 2020 amongst other studies (e.g.Salim & Lee 2012), when considering galaxies found at the bright-end of the rest-frame UV luminosity function (UVLF), the effect of scatter must be considered.If we consider an underlying stellar mass function (SMF) with an exponential cut-off above the characteristic mass, and estimate from this the expected UVLF, given the scatter between  ★ and the  UV , we expect to see a shallower decline in the number density to brighter galaxies.Because of the steepness of the SMF, galaxies that are found to be bright in the rest-frame UV are necessarily the sources that have low dust attenuation.Galaxies of the same stellar mass, but with a higher than average attenuation would instead be scattered fainter on the UVLF and be lost within the large population of fainter sources.The majority of the REBELS sources were selected from widearea ground-based data and they represent a very rare population of galaxies that sit brightward of the knee of the UVLF.Hence from the arguments detailed above, we would expect them to show a relatively low dust obscuration given their stellar mass (e.g. as found in Algera et al. 2023b), and hence a lower IRX and bluer .An additional factor that could contribute to the unexpectedly blue colours we observe for the REBELS sources is a geometric offset between the dust and stars within the galaxy (e.g. as predicted by Popping et al. 2017).The UV variability model of Shen et al. (2023) predicts that the rest-frame UV brightest galaxies should be blue, due to their young ages and clearance of dust in the star-burst phase.A clumpy morphology with an offset between the observed young stars stars and dust, leading to relatively unobscured regions, has already been observed in higherresolution ALMA observations of REBELS galaxies in Bowler et al. (2022) Figure 8.The modified blackbody fitting parameters from previous observations and models.We show the dust temperature against the emissivity coefficient, for an optically thin model.The results for the theoretical analysis by Sommovigo et al. (2022a,b) are shown as the purple circle (square) for REBELS (ALPINE).In the ALPINE survey, Béthermin et al. (2020) measured an empirical SED that had equivalent modified blackbody parameters as shown by the blue squares (the  = 4 point has the smaller error bar, with the  = 5 point showing a higher temperature).For both these studies the  d value was fixed.The contours show the offset computed in the derived IRX value from the parameters used in this work ( dust = 46 K,  d = 2.0).Such that the Béthermin et al. (2020) parameterisation would give a lower IRX by ∼ 0.25 dex.The results of fitting with a self-consistent (not optically thin) model to  = 4-8 galaxies with multiple FIR data points in Witstok et al. (2023) are shown as the grey open diamonds, demonstrating the current uncertainties and potential intrinsic scatter on these parameters at high redshift.

Implications for the dust attenuation curve at z > 4
We find that the IRX- relation derived from a stacking analysis of the ALPINE and REBELS samples are consistent with a Calzetti-like relationship as found at  = 0, with our assumed rest-frame FIR SED for both samples.We find no evidence that the sources on average lie below the relation, which could indicate a different attenuation law such as an SMC-like relation (discounting the effect of geometry which tends to make the observed attenuation law appear shallower even in the case of a steeper extinction law).As shown in Fig 7, our results differ from those found using the ALPINE dataset by Fudamoto et al. (2020b), who found a deficit in IRX for a given rest-frame UV slope and concluded that an SMC-like attenuation curve was preferred.As discussed in Section 4, the difference between our results and those of Fudamoto et al. (2020b) is primarily due to the assumed rest-frame FIR SED (with this work using a higher dust temperature by 5 K), with a minor effect of our analysis deriving bluer rest-frame UV slopes from the stacked photometry (Δ = 0.1).The best-fit relation presented in Fudamoto et al. (2020b) assumed an intrinsic  0 = −2.62 (fits were also presented for redder  0 ) and showed a tentative evolution in the normalisation, such that the IRX is decreasing at a given -slope from  = 4.5 to  = 5.5.We instead find little evolution in the IRX- relation from  ≃ 4 to  ≃ 7.2 from our analysis.This relies on our assumption that the FIR SED (and intrinsic rest-frame UV slope) remains approximately constant between the samples and redshift ranges, as found by Sommovigo et al. (2022a,b) which motivated the  d and emissivity coefficient used in this work for the ALPINE and REBELS datasets.If instead the dust temperature differs between the two samples at above and below  = 6, then we would conclude that the REBELS sources lie above the ALPINE galaxies on the IRX- plane by approximately 0.3 dex.This would potentially indicate a different selection function between the surveys (e.g.we already know that ALPINE targeted more Lyman- emitters than REBELS; see Section 5.3).
Our results are also higher (particularly at redder ) than the best-fit derived by Schouws et al. (2022) (assuming an intrinsic  0 = −2.23,at  d = 50 K).However the stacked results from that work are consistent within the errors with our findings from REBELS (see Fig. 5).In the range of stellar mass between log 10 ( ★ /M ⊙ )= 9-10, galaxies in both REBELS and ALPINE appear blue ( ≃ −2) and the obscured fraction is 0.4-0.6 (Table 2).For the higher mass galaxies found within ALPINE (as REBELS contains very few galaxies at log 10 ( ★ /M ⊙ )> 10 assuming a parametric SFH; although see Topping et al. 2022) we find redder slopes ( ≃ −1.5) and obscured fractions approaching 0.9.Comparing our results to studies that targeted lower mass galaxies there is some evidence for a difference in IRX with galaxy luminosity (or stellar mass).For example the stacking analysis of Bouwens et al. (2016) and Bouwens et al. (2020) derived from the deep ALMA data covering the Hubble Ultra Deep Field showed only upper limits over the broad redshift range  = 4-10.Bouwens et al. (2020) also found evidence for an SMC-like attenuation curve at log 10 ( ★ /M ⊙ )< 9.25, while a Calzetti like curve was preferred at higher masses.These results are consistent with our analysis, which supports a Calzetti-like attenuation curve for log 10 ( ★ /M ⊙ )> 9 from  = 4-8.Our find that there is a lack of evolution in the attenuation curve from local starburst galaxies up to  ≳ 7 could suggest that the conditions required to form Calzettilike dust are present already 800 Myr after the Big Bang.As shown by various theoretical studies however, taking the observed IRX- relation as evidence for a particular dust attenuation curve can be problematic due to the myriad of other factors that can influence the observed relation (e.g.Popping et al. 2017;Narayanan et al. 2018).The fact that we do not see a significant excess above the Calzetti relation further suggests that on average the REBELS and ALPINE galaxies are not dominated by significant optically thick regions contributing to the  IR , that would boost the measured IRX values (e.g. as seen in the case of ULIRGs; Casey et al. 2014).

Evidence for a high gas-phase metallicity?
Further insight into the origin of our results for massive galaxies can be gained from the works of Shivaei et al. (2020a,b) who found that the dust attenuation curve correlated most strongly with gasphase metallicity from a detailed study of sources at  = 2-2.5.They found that at 12 + log(O/H) > 8.5 a Calzetti-like attenuation curve was consistent with the data, with a steeper attenuation curve found for lower metallicity sources.Qualitatively our results agree with these  ≃ 2 studies, although the difference in methodology makes it non trivial to compare quantitatively.In particular, Shivaei et al. (2020b) remove sources with very strong emission line strengths ([OIII] 4959, 5007 > 630Å) that we inferred to be present in the majority of the REBELS sources (Bouwens et al. 2022).In addition, they compute the IR luminosity using a template from Rieke et al. (2009) for the more massive sources in their sample.These templates are stated to be equivalent to a greybody curve with  d = 38-64 K and  d = 0.7-1, which as shown in Fig. 8 would give lower values of the IRX by 0.3-0.4dex.With these caveats in mind, these previous results would suggest that the ALPINE and REBELS galaxies with log 10 ( ★ /M ⊙ )≃ 9.5 show an increased gas-phase metallicity compared to lower mass sources, that then leads to the observed Calzetti-like dust attenuation curve.If we take the mass-metallicity relation derived by Sanders et al. (2021) at  ≃ 3, we would predict a metallicity for this sample of 12 + log(O/H) ≃ 8.3.We also obtain a comparable estimate of the metallicity using the fundamental mass-metallicity relation (FMR) from Curti et al. (2020) assumming a total SFR ∼ 50 M ⊙ /yr (Table 2).Recent results from early JWST analysis for lower mass galaxies have revealed similar FMRs as at  ≃ 2 (Nakajima et al. 2023) as well as evidence for a deficit from the FMR in galaxies at  ≃ 7 (Curti et al. 2023).If this deficit is found to hold for galaxies at log 10 ( ★ /M ⊙ )≳ 9.5, this would imply a lower metallicity by ∼ 0.3 dex for galaxies in our sample given their  ★ .Regardless of which of these metallicity calibrations we consider, the derived values of average metallicity quantitatively disagree with the results of Shivaei et al. (2020b) that would predict that these sources should show an SMC-like dust attenuation curve at 12 + log(O/H) < 8.5.Despite this, our results are qualitatively in agreement with a picture where the ALPINE and REBELS sources have higher metallicities than lower-mass/fainter sources at the same redshifts, and hence a Calzetti-like dust attenuation law remains the preferred fit even up to  ≃ 7.With upcoming NIRSpec observations we will be able to directly determine the extension of the FMR relation at  = 7 to higher masses, and test whether the sources that show significant dust detections are metal rich in comparison to lower mass (log 10 ( ★ /M ⊙ )< 9) galaxies.

The IRX-M ★ relation at z = 4-8
In principle the IRX- ★ relation provides a more fundamental physical relationship than the IRX- plane, where the latter can be affected by large errors on the rest-frame UV slope measurement and geometric effects (e.g.see Faisst et al. 2017;Liang et al. 2019;Sommovigo et al. 2020;Ferrara et al. 2022).The stellar mass represents an integral of all past star-formation activity in the galaxy, and hence is expected at least qualitatively to be correlated with the production of dust (e.g. via stellar dust production and the overall metal enrichment of the ISM; Dayal et al. 2022).Dunlop et al. (2017) showed that stellar mass is a strong predictor of a FIR detection (and thus  IR ; see also McLure et al. 2018;Bouwens et al. 2020).Up to  ≃ 3 there has been in general a good agreement between previous studies of the IRX- ★ relation, with the results found to be approximately consistent to the local relation (e.g.McLure et al. 2018;Koprowski et al. 2018;Álvarez-Márquez et al. 2019; although see Reddy et al. 2018).At higher redshift, Fudamoto et al. (2020a) found a steeper slope of the IRX- ★ relation at  = 2.5-4 implying a reduced obscured SFR fraction and total SFR than the relations recovered at lower redshifts (at log 10 ( ★ /M ⊙ )< 11 when the relations cross), a trend that was also recovered in the ALPINE sample (Fudamoto et al. 2020b).In Fig. 9 we show the measured IRX- ★ points from our stacking analysis of the ALPINE and REBELS samples.We find an offset in the measurements from the  = 3 results of Koprowski et al. (2018) and Álvarez-Márquez et al. (2019).These studies use empirical templates in their derivation of the FIR luminosity, with the best-fitting templates used in Koprowski et al. (2018) showing a temperature of 40 K.Given the observed relation of  d with redshift, and the lack of information on the FIR SED at  > 6 and observed changes in the FIR SED with physical properties (e.g. ★ , SFR; Álvarez-Márquez et al. 2019), it is not trivial to compare IRX values over the redshift range  = 2-8 and be fully confident of the exact offsets between studies.Nevertheless, we recover an offset of between 0.5-1.0dex when comparing our  > 4 results to those at  ≃ 3, depending on the lower redshift relation assumed.Fitting these data points we derive the following IRX- ★ relation: log 10 (IRX) = 0.69(±0.12)log 10  ★ 10 10 M ⊙ + 0.40 (±0.05).(2) While the derived slope is only weakly constrained due to the small dynamic range in stellar mass that we probe with the ALPINE and REBELS samples, it is consistent with the  ≃ 3 relations.Further analysis of galaxies at log 10 ( ★ /M ⊙ )< 9 will be required to confirm if the slope does steepen at  ≳ 5. Due to the strong  ★ dependence of obscuration however, these measurements are extremely challenging and must rely on stacking (for example the ASPECS program only found 18 ALMA detections from a sample of 1362 galaxies at  = 1.5-10, with the highest redshift being at  = 3.7; Bouwens et al. 2020).Regardless of the exact form of the relation, the offset we observe in the derived IRX for a given  ★ shows that between  ≃ 2 and  = 4-8 the degree of obscured star formation has dropped by a factor of ∼ 3-10.This is despite the sources showing good agreement with the galaxy "main-sequence" (Topping et al. 2022;Algera et al. 2023b) and hence they do not show a deficit in total SFR.This finding is consistent with previous results that have derived a lower obscured fraction as a function of  ★ at  ≳ 5 (Fudamoto et al. 2020b;Gruppioni et al. 2020;Algera et al. 2023b), however we do not recover the steepness of these relations.We note that recent results on the dust obscuration as a function of stellar mass from Shapley et al. (2023) using the Balmer decrement measured with JWST have not found evidence for evolution in the degree of dust obscuration between  ≃ 2-6, in contrast to our results.One explanation for these differing conclusions is that there is an evolution in the relation between continuum (as we measure in this work using the rest-frame UV slope) and the nebular reddening (as measured by Shapley et al. 2023).Future rest-frame FIR observations of samples of massive  > 7 galaxies, such as those selected via upcoming wide-area NIR imaging observations from Euclid and Roman, will be required to confirm if the reduction in dust obscured star-formation continues to even higher redshifts.

Caveats
A key assumption in this work is that of a fixed dust SED with a dust temperature of  d = 46 K and  d = 2.0.These parameters were utilized within the REBELS survey as the fiducial rest-frame FIR SED and were derived by Sommovigo et al. (2022a).The dust temperature at  > 5 remains uncertain, and recent results have suggested that it may vary between galaxies from  d = 35-90 K (e.g.Hashimoto et al. 2019;Bakx et al. 2020;Algera et al. 2023a).For the emissivity index, Witstok et al. (2023) found a best-fitting value of  d = 1.8 ± 0.3, consistent with our assumed model.These results would suggest that a different dust temperature would be required for each source.However, given the large errors for individual sources, and the lack of these measurements for the full sample, this analysis is unfeasible at this time.
Another important caveat is that the ALPINE and REBELS samples were selected in different ways, and for example the  ≃ 5.5 sample from ALPINE includes a larger fraction of Lyman- emitters than REBELS.The REBELS survey has no confirmed Lyman- emission stronger than  0 = 25 Å (Endsley et al. 2022), while in the  = 5.5 bin of ALPINE 57 percent of the sources have  0 > 25Å.This comparison is complicated by the impact of the neutral IGM as we approach  ≃ 7, however as galaxies selected as Lyman- emitters have been shown to have lower dust content than LBGs (e.g.Schaerer et al. 2015;Matthee et al. 2017) this difference could cause a bias to lower IRX values in the ALPINE sample.There is also evidence at lower redshifts that the dust attenuation curve varies between galaxies of the same stellar mass (e.g.Salim & Boquien 2019; see the review by Salim & Narayanan 2020).Using a subset of the ALPINE sample, Boquien et al. (2022) derived a range of best-fitting attenuation curves also suggesting that the attenuation law varies significantly between individual galaxies.Narayanan et al. (2018) have predicted using a simulation that the variation in attenuation curves between galaxies decreases with increasing redshift, however further observations are needed to understand the properties of dust within high-redshift galaxies.
Another factor to consider, with reference to our inferences about stellar mass dependence, is that we assumed a certain parametric star-formation history in computing these parameters.Topping et al. (2022) has shown that these measurement may underestimate the mass by on average 0.5 dex.This would cause an increased deficit (up to around 0.35 dex, given the best-fit relation) in the IRX- ★ relation to previous studies, however care must be taken to compare similar methodologies in determining  ★ .Future observations of the REBELS galaxies with JWST (e.g. through program PID1626, PI Stefanon) will constrain the SFH and hence stellar masses with greater accuracy.
Finally, we consider the potential effects of sample incompleteness.Due to the rest-frame UV based Lyman-break selection (or Lyman- selection for some ALPINE sources), our samples will be incomplete to any highly obscured sources.Such sources have been identified with ALMA, for example the 'UV-dark' galaxies found serendipitiously via their [CII] emission only 40-60 pkpc from the central (rest-UV bright) REBELS source in Fudamoto et al. (2021).It is challenging to place these types of galaxies on the IRX- or IRX- ★ relation as their rest-frame UV and optical are undetected, however it is likely that they are extremely red and massive (log 10 ( ★ /M ⊙ )≲ 10-10.5)and with a high IRX.Further study is needed to understand this population of extremely obscured galaxies at  ≃ 7.

COMPARISON TO MODELS
The attenuation curve and the IRX- relation have been the subject of many theoretical and simulation analyses (e.g.see recent review by Salim & Narayanan 2020).Here we discuss the handful of results that focus specifically on predictions for the high-redshift Universe.In a work that was based on the REBELS sample, Dayal et al. (2022) provided the first theoretical comparison to the dust properties of the REBELS galaxies.Using the semi-analytic model (SAM) DELPHI, they were able to match the observed dust masses while also reconciling the intrinsic UV luminosity function with the observed (attenuated) luminosity function at  = 7. Ferrara et al. (2022) further identified that some REBELS sources show inconsistencies between the rest-frame UV properties () and FIR emission, suggestive of spatially offset regions within the galaxies.We recover this result in the form of sources that lie above the typical Calzetti-like IRX- relation -showing bluer than expected colours for their measured  IR .While Ferrara et al. (2022) argue that this makes the IRX- relation difficult to utilize for these galaxies, in this work we have found that on average there is a correlation in the IRX- plane for galaxies in the REBELS sample (that is consistent with a Calzettilike relation).In Fig. 10 we present the most recent analysis of dust emission from the DELPHI SAM in comparison to our derived IRX- results from REBELS and ALPINE.We took the output of the model presented in Mauerhofer & Dayal (2023) and computed the rest-frame UV slope using an identical method to that used on the ALPINE and REBELS stacks, by fitting a power law to the resulting predicted SED.Sources were extracted that had −24 <  UV < −20, and nebular continuum emission is included in the modelling of the SED.As is evident from the figure, comparing the DELPHI model of Mauerhofer & Dayal (2023) to the results of this work we find a good agreement with the model predictions.The predicted change in IRX and  with increasing stellar mass matches very well with what we find in the data, with the observed galaxies in ALPINE at log 10 ( ★ /M ⊙ )> 10 lying on top of the model predictions for this mass range.Furthermore, the bluest  values are comparable to those found in the data for galaxies at log 10 ( ★ /M ⊙ )≃ 9, and the lack of evolution we see between  = 4-8 is also recovered in DELPHI.As shown in Fig. 10, the SAM predicts that the IRX- should increase with redshift, due to the increased optical depth within compact galaxies at the highest redshifts.While there are assumptions made in the computation of the IRX- in both the models (e.g. through escape fraction of UV photons, dust heating etc.) and in the observations (with the systematic uncertainties in the  IR determination), it is reassuring that the general trends are present in this basic comparison.
Turning now to works that consider hydrodynamical simulations, Narayanan et al. (2018) used the MAFUSA model to predict the attenuation law at  = 6.They found a relatively grey attenuation curve, with the scatter between different galaxies becoming reduced at higher redshift due to more consistent stellar ages between sources and an increase in complexity in the star-dust geometry.The predicted curve is very close to that of the Calzetti et al. (2000) curve, however with a prominent 2175Å bump.Vĳayan et al. (2023) used the FLARES simulation to investigate the expected dust attenuation curve for galaxies at  = 5-10, with a particular focus on examining the effect of star-dust geometry.For galaxies with a similar SFR to  2).The results from DELPHI from  = 4.6 to  = 7.9 are shown as the small dots, with the colour scale corresponding to log 10 (  ★ /M ⊙ ) as derived in the model.There is a slight evolution within the redshift range, which we highlight with the solid and dotted grey lines from  = 7.9 and  = 4.6 respectively.
the REBELS and ALPINE sources, they find a dust law with a similar slope to that of the Calzetti et al. (2000) relation, despite the input of an SMC-like extinction law to their simulations.In FLARES, galaxies with a higher SFR show a larger degree of clumpiness, which in turn leads to a larger range in  V over the galaxy surface and hence a greyer attenuation curve when integrated across the source.The zoom-in cosmological SERRA simulations presented in Pallottini et al. (2022) produces galaxies at  ≃ 7.7 that are in excess of the majority of observations at high redshift (and in excess of the Calzettirelation).Pallottini et al. (2022) attribute the differences to first, the uncertain FIR SED and hence systemic errors in the derived  IR from observations (they suggest that a higher dust temperature should be assumed, exceeding 90 K), and second, to potential inaccuracies in the feedback prescription in the modelling and an insufficient spatial resolution to fully resolve the molecular clouds where the majority of the FIR luminosity is produced.Finally, Liang et al. (2019) used the cosmological simulation MassiveFIRE to model the IRX- relation of  = 2-6 galaxies.They found a relation similar to that found locally in the Calzetti-relation (although as a Milky Way dust attenuation curve was input, this is perhaps not surprising) Liang et al. (2019) directly predict the IRX- values for the galaxies in their simulations, finding objects that are blue ( < −1.5) and occupy the lower left region of the diagram.They conclude the scatter around the relation is dominated by different intrinsic  0 slopes.In conclusion, in general these models predict a dust attenuation law that is similar to the local Calzetti et al. (2000) relationship for galaxies at  > 6 with similar stellar masses and SFRs to the REBELS sample.When the models have predicted the IRX- relation this has been in agreement, or in excess of, what we observe for our combined ALPINE and REBELS samples (Liang et al. 2019;Pallottini et al. 2022;Mauerhofer & Dayal 2023).We therefore conclude that cur-rent models and observational constraints from this work favour a distribution in the IRX- plane that is consistent with predictions for a Calzetti-like dust attenuation law at  = 4-8, with no evidence for a deficit from the local relation towards what would be predicted by a screen of SMC-like dust.

CONCLUSIONS
In this study we present an analysis of the ALMA large program REBELS, which provided observations of the dust continuum emission for a sample of 49 galaxies at  = 6.5-8.We also perform a consistent analysis of the ALPINE large program that targeted galaxies in the redshift range of  = 4-6, to provide a key comparison to the REBELS results and investigate any evolution in the dust emission properties with redshift.The main conclusions of this work are as follows: • When compared to the expected colour from the extrapolated colour-magnitude relation for fainter galaxies, the REBELS sources are bluer than expected by up to Δ = 0.5-1.0(depending on the extrapolated relation).These results point to a flattening, or potential turnover, of the colour magnitude relation bright-ward of  UV = −21.5, that can be understand as a consequence of scatter on an underlying steep galaxy stellar mass function.In this scenario, the REBELS galaxies represent the sub-set of sources at a given  ★ that have low dust attenuation and hence have bright rest-frame UV magnitudes and a blue colour.
• When stacking the REBELS sources in the ALMA Band 6 data we find detections at > 5 significance for all but the highest redshift and brightest  UV bins.We derive the  IR and IRX for the sample from these results, assuming a modified blackbody curve with a dust temperature of 46 K and  d = 2.0 as derived in the model of Sommovigo et al. (2022a).The IRX- relation we derive for the REBELS sample with this assumed FIR SED is as expected for a Calzetti-like attenuation law (with an intrinsic rest-frame UV slope of  0 = −2.3).In comparison to other studies at  ≃ 7, we find no strong evidence for a systematic deviation below this relation, although large scatter is found between individual sources.
• By assuming the same FIR SED for stacks of the ALPINE data (as motivated by the study of Sommovigo et al. 2022b), we find that the IRX- results at  ≃ 4.5 and  ≃ 5.5 produce values that are consistent with our REBELS findings.We therefore find negligible evolution in the IRX- relation from  = 4-8, and conclude that there is little evidence for a deficit in the IRX- relation at  > 4 as compared to the local starburst results of e.g.Calzetti et al. (2000).Comparisons to previous studies at  > 4 again highlight that the assumed FIR SED has a dramatic affect on the derived  IR and care must be taken in comparing different studies.
• We compute the IRX- ★ relation from our combined ALPINE and REBELS samples, finding a similar slope to previous analyses at  = 3, but with a 0.5 dex offset to lower IRX at a given  ★ for our assumed FIR SED.These results corroborate previous findings (e.g.Schouws et al. 2022;Algera et al. 2023b) that show that for a given stellar mass, the proportion of obscured star-formation is reduced at  > 4 by a factor of ≳ 3.
• In comparison to models of dust attenuation in  ≳ 6 galaxies, we find that in general a Calzetti-like attenuation curve is predicted.In simulations that compute the IRX- relation we find that the results are in good agreement with our combined ALPINE and REBELS analysis, with few studies predicting very low 'SMC-like' relations.In a detailed comparison to the predictions of the DELPHI SAM, we find very good agreement between this model and our observations over the stellar mass range of log 10 ( ★ /M ⊙ )= 9-11.These results indicate that despite complexities and caveats in both the modelling and observation of high-redshift galaxies, qualitatively (and quantitively for the DELPHI model) there is good agreement between the predicted and measured trends of rest-frame UV colour, dust attenuated star-formation and stellar mass at  = 4-8.
To understand the origin of the effects seen in this work and others, a more detailed view on the measured properties of galaxies (such as more precise  ★ ,  and gas-phase metallicity) are required.Cycle 1 JWST observations of 12 of the REBELS sample are being obtained as part of the General Observer program PID1626 (PI Stefanon).This program will target the sources with the NIRSpec Integral Field Unit, providing gas-phase properties from the rest-frame optical emission lines for this unique sample.The data will also provide refined (and resolved) rest-frame UV slope and  ★ measurements.For the ALPINE survey, Cycle 2 observations have been approved for a subset of the most massive 18 galaxies for NIRSpec IFU from GO program 3045 (PI Faisst).The results of these programs, coupled with ongoing additional follow-up with ALMA to determine the dust temperature (via additional bands) and spatial distribution of the gas, dust and stars (via higher spatial resolution ALMA observations), will provide additional insights into the evolution of the IRX- and IRX- ★ relations at  = 4-8.

APPENDIX A: ALPINE ANALYSIS
We show the results of our analysis of the ALMA sample for the two redshift ranges centred on  ≃ 4.5 and  ≃ 5.5 in Fig. A1.The data are tabulated in Table A1.The individual points are plotted using the -values as derived by fitting to the best-fit SED by Faisst et al. (2020b) with some objects hitting the edge of parameter space, while the IRX value was recomputed in this work.The errors on the individual -slopes are smaller in ALPINE than in REBELS, due to the sources being in general brighter and having more photometric bands covering the rest-frame UV part of the SED.There are a larger fraction of non-detections in the dust continuum in ALPINE (80 percent; Fudamoto et al. 2020b) in comparison to REBELS (63 percent; Inami et al. 2022).From the bootstrap resampling used to The individual galaxy values where there is a detection in the dust continuum are shown as the filled grey points, with undetected sources shown as the light grey limits.The results of our stacking analysis, which takes into account the extended dust emission observed for the stacks at log 10 (  ★ /M ⊙ ) > 10, are shown as the blue and red points.The Calzetti and SMC-like relations are shown as described as in Fig. 4. derive errors on the stacked ALMA flux, we also find evidence for an increased scatter within the bins than found in the REBELS sample.
Table A1.The measured FIR fluxes and derived properties of the eight ALPINE stacks.The top (bottom) two rows show the lower (higher) redshift stack, with the rows ordered by  UV bin.Columns 1 and 2 detail the redshift and number of sources included in each stack.The average redshift and  UV of each stack are shown in Columns 3 and 4. Column 5 presents the measured peak ALMA flux, with the corresponding S/N shown in brackets.The flux measured using a Gaussian fit is shown in Column 6.The derived FIR luminosity (assuming a  dust = 46 K and  d = 2.0) and the resulting IRX values are shown in Columns 7 and 8.The  IR was determined from the peak flux for log 10 (  ★ /M ⊙ ) < 10 and the Gaussian flux at higher masses, as the log 10 (  ★ /M ⊙ )= 10-11 stack shows evidence for being spatially extended.Finally the rest-frame UV slope  is presented in Column 9, as measured from the best-fitting SED model.

Mass
Figure1.The average optical and NIR photometry and corresponding bestfitting SED model for the REBELS samples at  < 6.9 (top) and  > 6.9 (bottom).In each plot we show the observed fluxes and errors from the bright and faint stacks as the blue and red filled circles respectively.The six points correspond to the VISTA     and Spitzer/IRAC[3.6m]and  [4.5m]  bands.The open circles correspond to the synthetic photometry as Fig. 1.For ALPINE, we fit to the    [3.6m] [4.5m]bands for the  = 4.5 stack, and to the    [3.6m] [4.5m]bands for the  = 5.5 stack.A delayed- model was assumed (Φ() ∝  − (/ ) ) in which the timescale of the decline was allowed to vary in the range  = [0.3,10.0] Gyr and the age from 1 Myr up to the age of the Universe at that redshift.The metallicity was fixed to 0.2 Z ⊙ , and the Calzetti et al. (2000) dust law was assumed with the attenuation in the -band constrained to the range  V = [0, 2].We allowed the nebular ionization parameter to vary in the range log 10 () = [−2, −4].Uniform priors were assumed for all of the fitted parameters.These parameters resulted in acceptable fits to the REBELS sources as seen in Fig. 1, with no evidence for truncation of the resulting corner plots.Assuming a different SFH (e.g.constant or ) or metallicity only marginally affected the  ★ by at most 0.1 dex and the derived  values by < 0.05.Note that assuming a non-parametric SFH for the REBEL sample as presented in Topping et al. (

Figure 2 .
Figure 2. The ALMA Band 6 ( rest ≃ 150 m) stacks for REBELS.The stacks in the redshift range 6.5 <  < 6.9 (6.9 <  < 7.7) are shown in the upper (lower) row.The left and right columns show the brighter and fainter  UV stacks respectively, with the brighter stacks containing 7 (6) sources and the fainter stack containing 13 (14) sources at  ≃ 6.7 ( ≃ 7.2).The stamps are 10 arcsec on a side, with N to the top and E to the left.The colour scale is saturated beyond the range [ −2.0, 10.0]  and contours are shown at 1 intervals.The average beam for the data included in the stack is shown as the grey ellipse, with the position angle determined as the mode of the input image values.The stacks are consistent with being unresolved at the resolution of the data (1.2-1.6 arcsec FWHM).
Figure 3.The REBELS galaxies at  = 6.5-6.9 (upper) and 6.9 <  < 7.7 (lower) in comparison to the colour-magnitude relation found by previous studies.In each plot the individual galaxy measurements of the rest-frame UV slope and  UV are shown as the grey open points, while the stacked results (from the fits shown in Fig.1) are shown as the blue filled points.The relationship derived from fainter studies are shown as the black solid, dotted and dashed lines fromCullen et al. (2023),Topping et al. (2023) andBouwens et al. (2014) respectively.Slight differences between the relations in the upper and lower plot are due to the derived evolution of the relation from these studies.We also show the data points from the same three works (at a fixed rather than evolving redshift), note in particular that the redshift range of the data from the Cullen et al. (2023) study is at  > 8.

Figure 4 .
Figure4.The IRX- relation as derived from the REBELS sample.The individual dust continuum detected galaxies are shown as the filled dark grey points, with undetected galaxies shown as the lighter grey upper limits.The results from our stacking analysis are shown as the dark blue circles (light blue squares) at  ≃ 6.7 ( ≃ 7.2).Within each redshift bin the brighter  UV stack is found to have a redder  and a larger IRX.We assume a dust temperature of 46 K and emissivity of  d = 2.0, with the arrow shown in the lower right of the plot illustrating the systematic uncertainty we would obtain assuming ±15 K.The expected relation for Calzetti-like dust assuming  0 = −2.3 is shown as the black solid line, with the expected relation from an SMC extinction curve shown as the dashed line.The right axis shows the dust attenuation in magnitudes corresponding to the IRX according to the Calzetti dust law with a screen geometry.

Figure 5 .
Figure 5.The stacked IRX- points from our analysis of the REBELS sample in comparison to previous results for LBGs at  ≃ 7. We show the stacked results from Schouws et al. (2022) as the orange diamonds.The six individual galaxy results of Bowler et al. (2022) are shown as the open grey circles.Further sources from Witstok et al. (2022), Bakx et al. (2021) andMolyneux et al. (2022) are shown as the purple diamonds, red square and purple diamond respectively.We show the radio-loud AGN identified initially as a bright  ≃ 7 LBG byEndsley et al. (2023) as the black star in the upper right.

Figure 7 .
Figure 7.The IRX- relation derived in this study from the ALPINE and REBELS surveys at  ≃ 4-8.The ALPINE points at  = 4.5 ( = 5.5) are shown as the orange (red) open squares, while the REBELS points at  = 6.7 ( = 7.3) are shown as the blue (light blue) circles.Our best fitting IRX- relation to these data are shown as the solid blue (light blue) lines for an assumed intrinsic rest-frame UV slope of  0 = −2.3(−2.5).We also show the best-fitting IRX- relation from Schouws et al. (2022) (assuming a  0 = −2.23)as the blue dotted line.The best-fitting IRX- relations assuming a  0 = −2.62 from Fudamoto et al. (2020b) are shown as the orange (red) dashed lines for  = 4.5 ( = 5.5).The difference between our ALPINE results and those of Fudamoto et al. (2020b) can be mostly explained by the different assumed FIR SED.Fudamoto et al. (2020b) assumed  d = 41 K, in comparison to the  d = 46 K assumed in this work, leading to the difference illustrated in the bottom right as the red arrow.

Figure 9 .
Figure 9.The derived IRX- ★ values from our consistent analysis of the REBELS and ALPINE datasets.The results from our REBELS analysis are shown as the dark blue (light blue) circles at  ≃ 6.7 ( ≃ 7.2).The ALPINE results are shown as orange (red) open squares for the  ≃ 4.5 ( ≃ 5.5) sub-samples.The best-fitting relation to these data is shown as the black solid line with the grey shading indicating the 1 confidence regime.Relations from previous studies at  = 3 from Koprowski et al. (2018) and Álvarez-Márquez et al. (2019) are shown as the dotted blue and solid light blue lines respectively.The  ≃ 4-6 relation found by Fudamoto et al. (2020b) is shown as the dashed red line.

Figure 10 .
Figure 10.The predictions of the DELPHI semi-analytic model of Mauerhofer & Dayal (2023) for the IRX- relation, in comparison to our results from REBELS and ALPINE.The data points and lines are as shown in Fig. A1.The observed points from ALPINE that are found at  > −1.75 correspond to the log 10 (  ★ /M ⊙ )= 10-11 stacks of these data.The REBELS sources show log 10 (  ★ /M ⊙ )≃ 9.5 (Table2).The results from DELPHI from  = 4.6 to  = 7.9 are shown as the small dots, with the colour scale corresponding to log 10 (  ★ /M ⊙ ) as derived in the model.There is a slight evolution within the redshift range, which we highlight with the solid and dotted grey lines from  = 7.9 and  = 4.6 respectively.

Figure A1 .
Figure A1.The results of our reanalysis of the ALPINE dataset on the IRX- relation, with the upper (lower) plot showing the  ≃ 4.5 ( ≃ 5.5) results.The individual galaxy values where there is a detection in the dust continuum are shown as the filled grey points, with undetected sources shown as the light grey limits.The results of our stacking analysis, which takes into account the extended dust emission observed for the stacks at log 10 (  ★ /M ⊙ ) > 10, are shown as the blue and red points.The Calzetti and SMC-like relations are shown as described as in Fig.4.

Table 2 .
Chabrier (2003)023b2014)ividually results for REBELS presented inBouwens et al. (2022), the galaxies are massive, with  ★ ≳ 10 9 M ⊙ , and moderate ages of the order of 40-130 Myrs.The derived dust attenuation is relatively low, as expected from the fact that the sample is rest-frame UV selected and shows blue  slopes (see Section 4.2).Our most significantly FIR detected stack has the reddest  = −1.8 and strongest  V = 0.65 ± 0.05.We measure the unobscured SFR from the rest-frame UV emission using the luminosity at 1500Å and the conversion ofMadau & Dickinson (2014)for a constant SFR in the previous 100 Myr and a fixed metallicity of  = 0.1 Z ⊙ .The SFR from the FIR was derived using the conversion based on the same assumptions on the SFH fromMadau & Dickinson (2014), which provides an identical calibration to that used in the previous REBELS work byAlgera et al. 2023b.Both calibrations were adjusted to aChabrier (2003)initial mass function (IMF).The total SFR measured as the sum of these two components (SFR UV + SFR IR ) is in good agreement with that derived from the SED fitting.
H + [OIII] 4959, 5007 lines sit within a single [3.6m] or [4.5m] band.As shown in Fig. 1 we find that the REBELS sources are blue in the rest-frame UV, as probed by the    bands, with strong [3.6m]-[4.5m]colours evident.As we move to  > 6.9 we see a change in the IRAC colour indicative of H+[OIII] 4959, 5007 moving into the [4.5m]-band.The derived SED fitting parameters from this photometry including stellar mass and age are presented in supporting this hypothesis.