We study 31 little red dots LRD detected by JADESNIRCam and covered by the SMILESMIRI survey of which 70 are detected in the two bluest MIRI bands and 40 in redder MIRI filters. The medianquartiles redshifts are inlineformula mmlmath overflowscrollmmlmizmmlmimmlmommlmommlmsubsupmmlmrowmmlmn6.9mmlmnmmlmrowmmlmrowmmlmn5.9mmlmnmmlmrowmmlmrowmmlmn7.7mmlmnmmlmrowmmlmsubsupmmlmath inlineformula 55 spectroscopic. The spectral slopes flatten in the restframe nearinfrared consistent with a 1.6 m stellar bump but bluer than direct pure emission from active galactic nuclei AGN tori. The apparent dominance of stellar emission at these wavelengths for many LRDs expedites stellar mass estimation the medianquartiles are inlineformula mmlmath overflowscrollmmlmilogmmlmimmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmsubmmlmrowmmlmo stretchytruemmlmommlmrowmmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmsubmmlmommlmommlmsubsupmmlmrowmmlmn9.4mmlmnmmlmrowmmlmrowmmlmn9.1mmlmnmmlmrowmmlmrowmmlmn9.7mmlmnmmlmrowmmlmsubsupmmlmath inlineformula. The number density of LRDs is 10SUP4.00.1SUP MpcSUP3SUP accounting for 14  3 of the global population of galaxies with similar redshifts and masses. The restframe nearmidinfrared 24 m spectral slope reveals significant amounts of warm dust bolometric attenuation 34 mag. Our spectral energy distribution modeling implies the presence of lt0.4 kpc diameter knots heated by either dustenshrouded OB stars or an AGN producing a similar radiation field obscured by AV gt 10 mag. We find a wide variety in the nature of LRDs. However the bestfitting models for many of them correspond to extremely intense and compact starburst galaxies with massweighted ages 510 Myr very efficient in producing dust with their global energy output dominated by the direct in the flat restframe ultraviolet and optical spectral range and dustrecycled emission from OB stars with some contribution from an obscured AGN in the infrared.
['Galaxy formation', 'Galaxy evolution', 'High-redshift galaxies', 'Galaxy stellar content', 'Stellar populations', 'Broad band photometry', 'Galaxy ages', 'James Webb Space Telescope', 'Active galactic nuclei', '595', '594', '734', '621', '1622', '184', '576', '2291', '16', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Instrumentation and Methods for Astrophysics']
What Is the Nature of Little Red Dots and what Is Not MIRI SMILES Edition
{'What is the nature of Little Red Dots and what is not, MIRI SMILES edition': "Pablo G. P'erez-Gonz'alez , 1 Guillermo Barro , 2 George H. Rieke , 3 Jianwei Lyu , 3 Marcia Rieke , 3 Stacey Alberts , 3 Christina C. Williams , 4 Kevin Hainline , 3 Fengwu Sun , 3 D'avid Pusk'as , 5, 6 Marianna Annunziatella , 1 William M. Baker , 5, 6 Andrew J. Bunker , 7 Eiichi Egami , 3 Zhiyuan Ji , 3 Benjamin D. Johnson , 8 Brant Robertson , 9 Bruno Rodr'ıguez Del Pino , 1 Wiphu Rujopakarn , 10, 11 Irene Shivaei , 1 Sandro Tacchella , 5, 6 Christopher N. A. Willmer , 3 and Chris Willott 12  \n1 Centro de Astrobiolog'ıa (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrej'on de Ardoz, E-28850, Madrid, Spain 2 Department of Physics, University of the Pacific, Stockton, CA 90340 USA  \n3 Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA 4 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA 5 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK 6 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK 7 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK 8 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA 9 Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA 10 National Astronomical Research Institute of Thailand, Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand  \n11 Department of Physics, Faculty of Science, Chulalongkorn University, 254 Phayathai Road, Pathumwan, Bangkok 10330, Thailand 12 NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada", 'ABSTRACT': 'We study 31 little red dots (LRD) detected by JADES/NIRCam and covered by the SMILES/MIRI survey, of which ∼ 70% are detected in the two bluest MIRI bands and 40% in redder MIRI filters. The median/quartiles redshifts are z = 6 . 9 7 . 7 5 . 9 (55% spectroscopic). The spectral slopes flatten in the restframe near-infrared, consistent with a 1.6 µ m stellar bump but bluer than direct pure emission from active galactic nuclei (AGN) tori. The apparent dominance of stellar emission at these wavelengths for many LRDs expedites stellar mass estimation: the median/quartiles are log M ⋆ / M ⊙ = 9 . 4 9 . 7 9 . 1 . The number density of LRDs is 10 -4 . 0 ± 0 . 1 Mpc -3 , accounting for 14 ± 3% of the global population of galaxies with similar redshifts and masses. The rest-frame near/mid-infrared (2-4 µ m) spectral slope reveals significant amounts of warm dust (bolometric attenuation ∼ 3 -4 mag). Our spectral energy distribution modeling implies the presence of < 0 . 4 kpc diameter knots, heated by either dustenshrouded OB stars or an AGN producing a similar radiation field, obscured by A(V) > 10 mag. We find a wide variety in the nature of LRDs. However, the best-fitting models for many of them correspond to extremely intense and compact starburst galaxies with mass-weighted ages 5-10 Myr, very efficient in producing dust, with their global energy output dominated by the direct (in the flat rest-frame ultraviolet and optical spectral range) and dust-recycled emission from OB stars with some contribution from an obscured AGN (in the infrared).  \nKeywords: Galaxy formation (595) - Galaxy evolution (594) - High-redshift galaxies (734) - Stellar populations (1622) - Broad band photometry (184) - Galaxy ages (576) - JWST (2291) - Active galactic nuclei(16)', '1. INTRODUCTION': "In the very first month of JWST science operations, Labb'e et al. (2023a) identified a sample of compact sources with distinct spectral energy distributions (SEDs) presenting two clear breaks (Lyman and Balmer) in the NIRCam+HST data. The SEDs were also char-  \ncterized by a change of slope: they had red colors at observed wavelengths between ∼ 2 and (at least) ∼ 5 µ m, the range covered by the NIRCam long-wavelength (LW) channels, and a flat SED in the short wavelength (SW) bands. Labb'e et al. (2023a) identified these SWblue+LW-red sources as very massive, M ∗ > 10 10 M ⊙ , maybe significantly-obscured, A(V) > 1 . 5 mag, galaxies  \nat z = 7 -9. They have been claimed to represent a new galaxy population revealed through the new capabilities provided by JWST, easily reaching magnitudes down to 28-29 mag from 1 to 5 µ m.  \nThe existence of such early massive galaxies was quickly identified as a possible severe problem in our understanding of how galaxies form and evolve. Indeed, the large masses would be very difficult to explain with the current ΛCDM paradigm, since the amount of baryons already transformed into stars could exceed their abundance in early halos (Haslbauer et al. 2022; Boylan-Kolchin 2023; Dekel et al. 2023).  \nThis high-mass interpretation was challenged by adopting models with prominent nebular emission, which could account for their red nature due to the non-negligible contribution of emission lines (and nebular continuum) to the broad-band fluxes as they enter the NIRCam passbands for different redshifts, implying significantly smaller (by a factor of 10 or more) stellar masses (see Endsley et al. 2023a,b; P'erez-Gonz'alez et al. 2023a; Desprez et al. 2023). The estimated masses of early galaxies would also be reduced significantly if they had top-heavy stellar initial mass functions (Woodrum et al. 2023; Wang et al. 2023).  \nA simpler explanation in some cases is redshift errors or mis-identifications. Indeed, Kocevski et al. (2023) presented spectroscopy for one of the claimed massive galaxies in Labb'e et al. (2023a), showing an overestimation of the photometric redshift which could partly explain the high mass value for this and some other sources, as opposed to the more accurate photometric redshift and (10 times) smaller mass found in P'erezGonz'alez et al. (2023a) for this same source. However, the general agreement among different authors in the determination of photometric redshifts for this type of source implies that redshift errors are not worrisome in a statistical sense (estimations agree well between papers for ∼ 75% of the samples). In addition, some purported LRDs (at the 10-20% level) have been found to be brown dwarfs based on NIRCam colors, proper motions, and spectroscopy (Hainline et al. 2023a; Langeroodi & Hjorth 2023; Burgasser et al. 2023).  \nJointly with their colors, the extreme compactness of the bona fide high-redshift objects led to the term Little Red Dots (LRDs, as first suggested by Matthee et al. 2023b). Understanding LRDs has become more complex given that there is a variety of selection techniques that arrive at very similar types of objects in terms of compact morphology and red LW-to-SW colors, with high levels of overlap but also 'contamination' from other types of sources (e.g., not so compact red galaxies or little not-so-red dots). Labb'e et al. (2023a) originally  \nlooked for sources with SED breaks, which resulted in a flat or blue SED at wavelengths shorter than ∼ 2 µ mand a very red one at longer wavelengths to ≥ 5 µ m (SWblue+LW-red SEDs). Eventually arriving at similar objects, Barro et al. (2023), Akins et al. (2023), or Labb'e et al. (2023b) used colors such as F277W-F444W in their selection, adding additional constraints in the SW channels. Some LRDs also entered in the selection carried out by P'erez-Gonz'alez et al. (2023a), Barrufet et al. (2023) and Williams et al. (2023a) based on F150WF356W and F150W-F444W (see also red sources selected with F200W-F444W in Rodighiero et al. 2023 and Bisigello et al. 2023), which were in some cases included in larger investigations about the nature and cosmic relevance of 'HST-dark' sources. Finally, complementary to purely photometric techniques, spectroscopic data with less comprehensive photometry have also been able to identify LRDs, mainly based on the selection of strong H α or [OIII] emitters in blind NIRCam/grism spectroscopic surveys (Matthee et al. 2023a,b).  \nThe question about the general nature of LRDs was further tangled by finding evidence for the presence of Active Galactic Nuclei (AGN). Kocevski et al. (2023) identified a broad (1000 km s -1 ) component in the H α emission of one galaxy in the Labb'e et al. (2023a) sample, which implied the presence of a 10 7 M ⊙ supermassive black hole (SMBH). This large SMBH mass is consistent with the discoveries of high-redshift AGN in other studies, all indicating that black holes grew rapidly in the early Universe (D'ıaz-Santos et al. 2021; Larson et al. 2023; Maiolino et al. 2023a; Matthee et al. 2023b; Ubler et al. 2023; Harikane et al. 2023; Natarajan et al. 2023; Bogd'an et al. 2023), many lying above the Magorrian et al. (1998) relationship (Pacucci et al. 2023; Stone et al. 2023). Further spectroscopic analyses identified AGN in some additional LRDs and stellar emission in others, even both at different wavelengths in the same LRD (Matthee et al. 2023b; Greene et al. 2023; Killi et al. 2023). ALMA non-detections have also been used to defend the AGN nature of the bulk of the LRD emission (Labb'e et al. 2023b).  \nSolving the challenges outlined above requires going further into the mid-infrared where stellar emission and AGNcan be distinguished. MIRI data at observed wavelengths shorter than 8 µ m for a handful of LRDs was presented in Akins et al. (2023), Barro et al. (2023, also with a F1000W detection), indicating that they could be either dusty starbursts (see also Rodighiero et al. 2023) or obscured AGN, but no definitive conclusions could be reached with just the bluest MIRI bands. Williams et al. (2023a) used longer wavelength MIRI data from the SMILES program and found that the averaged SEDs  \n(stacked in observed bands) of the LRDs flattened beyond 5 µ m, interpreting this behavior as the expected turnover of a normal stellar SED at rest ∼ 1 . 6 µ m. In addition, and at odds with the AGN interpretation of the ALMA data of LRDs in the A2744 field presented by Labb'e et al. 2023b, Williams et al. (2023a) found that even deeper ALMA observations of LRDs in GOODSS could agree with SED models mostly dominated by stars (with some implications for the global dust emission SED shape).  \nAs is clear from the preceding discussion, the exact nature of the LRDs is still unclear, and there is no definitive proof that they are a single type of phenomenon, i.e., they could be a mixture of the various hypotheses advanced, their nature being also affected by selection biases.  \nIn this paper, we investigate the nature of LRDs using the best multi-wavelength, large area MIRI data available to date, those coming from the SMILES survey (Lyu et al. 2023). MIRI, in fact, is essential to disentangle the nature of LRDs. If dominated by star formation, the SEDs of LRDs should show a 1.6 µ m bump, redshifted to ∼ 10 -15 µ m at z = 5 -8 (e.g. Williams et al. 2023a). On the other hand, an obscured AGN will yield a steep SED. The significant uncertainty in the relative importance of star formation and nuclear activity needs both spectroscopic (e.g., Kocevski et al. 2023; Killi et al. 2023) and photometric data to cover the widest spectral range possible, as well as a careful SED analysis of a statistically representative sample. In modeling these sources, we will also explore the uncertainties by using four different sets of modeling software with a variety of approaches and assumptions.  \nThe main goals of this paper are to select and study LRDs in a comprehensive way, including examples representing a range of properties, and to see how our conclusions about them are influenced by different elaborated modeling approaches. To achieve these ends, the paper is organized as follows. Section 2 presents the NIRCam and MIRI data used to select and characterize a sample of LRDs in the GOODS-S field. Section 3 describes the average SEDs of LRDs and Section 4 presents stellar and AGN emission models that will be used to characterize the physical properties of LRDs on galaxy-by-galaxy and statistical bases in Section 5. Finally, Section 6 discusses and summarizes our conclusions. We also include in this paper three appendices that describe in detail the MIRI data reduction, the dust emission models that are the most important ingredients of the SED modeling, and provide SED fits for the whole sample.  \nThroughout the paper, we assume a flat cosmology with Ω M = 0 . 3 , Ω Λ = 0 . 7, and a Hubble constant  \nH 0 = 70kms -1 Mpc -1 . We use AB magnitudes (Oke & Gunn 1983). All stellar mass and SFR estimations assume a universal Chabrier (2003) IMF, which is a very relevant assumption for these parameters (but might not be accurate given the extreme conditions).", '2.1. NIRCam-based selection': "The sample of LRDs analyzed in this paper has been gathered from the JWST Advanced Deep Extragalactic Survey, JADES (Eisenstein et al. 2023a) Data Release 2 catalogs (Eisenstein et al. 2023b), which also gathered data from the First Reionization Epoch Spectroscopically Complete Observations (FRESCO, Oesch et al. 2023) and the JWST Extragalactic Mediumband Survey (JEMS, Williams et al. 2023b). We also used the Rieke et al. (2023) catalog for spectroscopic redshifts, which included values measured with data taken by NIRCam (Oesch et al. 2023) and NIRSpec (Bunker et al. 2023). Given the small (nearly pointlike) nature of our sources of interest and the fact that eventually we wanted to analyze SEDs including MIRI data, we used 0.25 '' radius PSF-matched photometry as our fiducial aperture for selection, but also checked the results using smaller radii, more suitable for NIRCam only. We searched for NIRCam SWblue+LW-red sources defined by F277W-F444W > 1 mag and F150W-F200W < 0.5 mag colors, and magnitudes F444W ≤ 28 mag, following the strategies presented in Barro et al. (2023), Labb'e et al. (2023b) and Greene et al. (2023). The two latter references further use a concentration criterion to select compact sources, but we found no need for it after applying the F150WF200W < 0.5 mag cut. In any case, the typical concentration of our original sample is indicated by the flux ratio (after applying the appropriate aperture corrections) of F(F444W) r =0 . 25 '' /F(F444W) r =0 . 10 '' = 1 . 04 1 . 10 0 . 98 (median and quartiles). That is, most of the sources are nearly point-like.  \nThe selection procedure is summarized in Figure 1, where we compare our sample with others found in the literature. Overall, our color selection criterion is bluer than that used in Barro et al. (2023), who imposed a redder color cut (F277W-F444W > 1.5 mag instead of 1.0 mag), and very similar to the one used for the sample presented in Labb'e et al. (2023b), although we impose a magnitude cut ∼ 0 . 5 mag deeper in F444W.  \nWe further restrict the sample to the area covered by MIRI within the Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES, Lyu et al. 2023), arriving at a sample of 37 galaxies. As shown by Hainline et al. (2023a), this type of color selec-  \nFigure 1. The left and central panels show the F 277 W -F 444 W vs . F 444 W color-magnitude and F 277 W -F 444 W vs. F 150 W -F 200 W color-color diagrams, as well as histograms of NIRCam colors, indicating the selection thresholds for LRDs ( F 277 W -F 444 W > 1 mag, F 150 W -F 200 W < 0 . 5 mag and F 444 W < 28 mag; dashed lines) relative to the bulk of the JADES DR2 galaxy catalog. The different colors indicate the subsets of LRDs detected in different MIRI bands: up to F1280W and beyond (Golden Five galaxies), up to F1000W, up to F770W, or not detected in MIRI (at the SMILES depth). Comparison LRD samples from Barro et al. (2023) and Labb'e et al. (2023b) in the CEERS and UNCOVER fields are shown with squares and triangles, respectively. The right panel shows the color vs. redshift diagram for the LRDs and JADES galaxies and the redshift distribution of LRDs (including median and quartiles). The 55% of the LRDs in our paper that have secure NIRSpecand NIRCam-based spectroscopic redshifts are marked with a black dot.  \n<!-- image -->  \ntion of high redshift galaxies might be contaminated by brown dwarfs. To account for that, and following their analysis, we removed from the original sample the 4 sources with color F115W-F150W < -0.5 mag. Two other LRDs were removed since they were identified with brown dwarfs with proper motion (Hainline et al. 2023a). Our final sample contains 31 galaxies detected in the SMILES 34.9 arcmin 2 area, i.e., the surface density is 0.9 LRD arcmin -2 . Postage stamps in RGB format constructed with NIRCam data are shown in Figure 2. Table 1 provides all relevant information about the selection of our sample of LRDs.  \nOur final sample includes 14 sources studied in Williams et al. (2023a), as well as 2 sources included in the sample of high redshift obscured AGN candidates presented in Lyu et al. (2023), one source identified as an AGN in Matthee et al. (2023b), and 5 galaxies selected as high-redshift candidates in Hainline et al. (2023b). Compared to the galaxies in common with Williams et al. (2023a), the sample in this paper is 0.5 mag fainter in F444W (medians and quartiles 26 . 5 27 . 5 26 . 0 mag vs. 26 . 0 26 . 9 25 . 4 mag), and bluer in F277WF444W (+1 . 3 +1 . 7 +1 . 1 mag vs. +1 . 6 +2 . 5 +1 . 5 mag) and F150WF200W (+0 . 0 +0 . 2 -0 . 1 mag vs. +0 . 2 +0 . 3 +0 . 0 mag). The sources in common with Hainline et al. (2023b) are among the faintest (median F444W 27.1 mag), slightly bluer (me-  \ndian F277W-F444W and F150W-F200W 1.1 mag and 0.0 mag, respectively), and at the highest redshifts in our sample (all with z > 8 . 2, median 8.4).", '2.2. Mid-infrared properties': 'In this section, we discuss the detections of our sample of LRDs in the MIRI data gathered by SMILES (Lyu et al. 2023). We first describe briefly this dataset, and then discuss the detection fractions in MIRI as a function of wavelength.', '2.2.1. MIRI data description': "The MIRI data gathered by the SMILES survey (program ID 1207, PI: G. Rieke, Lyu et al. 2023) were reduced with JWST pipeline version 1.12.3, reference files in pmap version 1138, which includes improved absolute photometric calibration and aperture corrections based on a better characterization of the PSFs, especially at the shortest wavelengths.  \nApart from the official pipeline steps, our reduction includes a special treatment of the background to deal with artifacts seen in the MIRI data (e.g., striping, tree rings, and inhomogeneities in rows and columns). This is achieved with a super-background strategy that was explained briefly in ' Alvarez-M'arquez et al. (2023) and Lyu et al. (2023). Given the importance of reaching the deepest fluxes possible in the analysis of high redshift  \n<!-- image -->  \nFigure 2. NIRCam F150W+F277W+F444W color composite postages of the sample of 31 LRDs in this paper. North is up, East is left, sizes are 2.5 × 2.5 arcsec 2 . Galaxies detected by MIRI at wavelengths longer than F1280W (Golden Five galaxies), F1000W, and F770W are marked in gold, red, and purple (and the sources are ordered according to this MIRI detections). Those not detected in MIRI are marked in magenta. We display redshifts, with spectroscopic values written with 4 decimals.  \nsources presented in this paper, we explain and characterize in detail the bespoke MIRI reduction procedures in Appendix A.  \nPhotometry in the MIRI bands was performed in circular apertures with radii of 0.2 '' , 0.3 '' , 0.4 '' , and 0.5 '' , applying the corresponding aperture corrections to each measurement. Uncertainties were calculated following the procedure explained in P'erez-Gonz'alez et al. (2023b), which gathers non-contiguous (i.e., independent) pixels to avoid the effect of correlated noise. Different measurements for our sample of LRDs agreed within 0.1 mag and we eventually chose for each band the weighted mean of all fluxes, typically the one corresponding to ∼ 70% encircled energy for each band. The photometry was revised visually to avoid contamination from cosmic ray showers, which we found to affect one source. We only kept fluxes with S/N > 5 and considered 5 σ upper limits otherwise.", '2.2.2. MIRI detections': 'Out of the 31 LRDs selected with the JADES NIRCam data down to F444W=28 mag, we detect [19,22,12,7,4,2,1,0] galaxies for MIRI bands [F560W,F770W,F1000W,F1280W,F1500W,F1800W, F2100W,F2550W] respectively. The corresponding detection fractions (in %) are [61,71,39,23,13,7,3,0] at > 5 σ levels of [26.1,25.8,24.7,24.3,24.2,23.0,22.6,20.8] mag (these magnitudes are measured in circular apertures of radius equal to the PSF FWHM, corrected to infinite extent). The detection fractions depend strongly on F444W magnitude. For the brightest galaxies, F444W < 26.5 mag, the detection fractions are (in %) [93,100,73,47,20,13,7,0] in [F560W,F770W,F1000W,F1280W,F1500W,F1800W, F2100W,F2550W], respectively.  \nWe divided the sample of LRDs into 4 subsets detected up to a given MIRI band: (1) 5 galaxies detected in F1280W and longer wavelengths, hereafter the Golden Five; (2) 7 galaxies detected only up to F1000W (i.e., excluding the previous type); (3) 7 sources detected only up to F770W (excluding the previous types); (4) and, lastly, the 12 galaxies not detected in any MIRI band. The panels and histograms in Figure 1, as well as the postage stamps in Figure 2, show these subsets with different colors. The first figure illustrates that the MIRI detection fraction depends strongly on the F444W magnitude and it is nearly independent of the F150WF200W color. Interestingly, the F277W-F444W color does exhibit some differences between subsets, with the Golden Five being among the bluest, and the F1000W sample among the reddest (we will comment on this in the following sections). It is worth mentioning that only  \ntwo of the Golden Five galaxies would be identified as LRDs under a more restrictive color selection of F277WF444W > 1.5 mag.', '2.3. Redshifts': "Atotal of 18 galaxies out of the 31 in the whole sample (55%) have spectroscopic redshifts provided by CANDELS (Guo et al. 2013), FRESCO (Oesch et al. 2023) data (measured by the JADES team) and JADES NIRSpec data (Bunker et al. 2023; Eisenstein et al. 2023b), all of them flagged as highly secure.  \nPhotometric redshifts for all galaxies were estimated with eazy (Brammer et al. 2008), using either direct flux measurements for all bands or with a modified version of the program that penalizes solutions implying fluxes above upper limits. We used v1.3 templates, which include a dusty galaxy and a high-EW emissionline galaxy spectrum, and we added a new extreme emission-line galaxy template as well as an AGN+torus model (based on the spectrum in Killi et al. 2023). Based on the comparison with spectroscopic values, we chose the maximum χ 2 photometric redshift as our main solution. We estimated redshifts both using only NIRCam fluxes and also adding MIRI measurements, obtaining better results when using NIRCam only (a conclusion based on the comparison with spectroscopic values).  \nThe quality of the photometric redshifts, when compared to spectroscopic values, is remarkably high, partly as a consequence of the strong emission lines present in many of the LRDs in our sample. The average absolute difference between photometric and spectroscopic values is 0.01, and we do not have any outliers (defined as galaxies with ∆( z ) / (1 + z ) > 0 . 1).  \nA histogram of the redshift distribution of our sample is given in Figure 1. The distribution is relatively flat between z ∼ 5 and z ∼ 9, peaking at z ∼ 7. This behavior is probably related to the fact that strong emission lines entering the F444W filter affect the selection significantly (as also seen in extreme emission line galaxies included in P'erez-Gonz'alez et al. 2023a). These lines are mainly [OIII]+H β and H α +[NII]+[SII], which translate to z ∼ 8 and z ∼ 6 for the lines lying within the F444W passband. Summarizing, the typical redshift (median and quartiles) of LRDs selected down to F444W ∼ 28 mag is z = 6 . 9 7 . 7 5 . 9 .", '2.4. Additional SED constraints': 'Apart from the NIRCam and MIRI fluxes, we used Hubble ACS fluxes measured in Hubble Legacy Field images (Illingworth et al. 2016; Whitaker et al. 2019). Except for one source, not covered by the NIRCam SW  \nchannels, we did not use WFC3 data. We checked that none of the LRDs in this paper are detected by Spitzer with the MIPS instrument or by Herschel with PACS and SPIRE. Taking into consideration the catalogs used in Barro et al. (2019), we imposed in our SED analysis 5 σ upper limits of 30 µ Jy for the MIPS 24 µ m band, and 4, 5, 9, 11, and 11 mJy for Herschel bands at 100, 160, 250, 350, and 500 µ m, respectively.', '3.1. Average spectral energy distributions': "Figure 3 shows the rest-frame SEDs, normalized at 0.4 µ m, and their average (more specifically, the weighted mean every 10 points calculated on a wavelength-ordered flux list) for the LRDs detected by MIRI and those with just upper limits. The data are compared qualitatively with some relevant templates to understand what the MIRI observations are telling us about the general nature of LRDs.  \nOverall, the plots highlight very clearly the characteristic flat, blue UV continuum plus steep, red optical continuum of LRDs. Such peculiar SEDs are difficult to model with traditional stellar population models under conventional assumptions and therefore require more complex or unusual treatment (to avoid biases in the determination of stellar masses, for example). Recent works have put forward different solutions to try to explain the emission of the LRDs using: strong, high EW( > 1000 ˚ A) emission lines that boost even the broadband photometry (Endsley et al. 2023a,b; Matthee et al. 2023b; Furtak et al. 2023; P'erez-Gonz'alez et al. 2023a; Desprez et al. 2023), complex stellar models with a flexible treatment of the dust attenuation (Labb'e et al. 2023a; Barro et al. 2023; Akins et al. 2023; Williams et al. 2023a), or hybrid models that combine stellar and AGN-driven emission with either component dominating different spectral regions or the full SED (Kocevski et al. 2023; Barro et al. 2023; Labb'e et al. 2023b; Greene et al. 2023).  \nThe fully AGN-dominated models have received more attention lately based on the spectroscopic observations of LRDs showing strong, and sometimes broad emission lines (FWHM > 2000 km/s; Kocevski et al. 2023; Harikane et al. 2023; Matthee et al. 2023b; Greene et al. 2023; Maiolino et al. 2023b). Briefly, these models combine AGN continuum emission from 1) a highly obscured accretion disk, 2) a small fraction ( ∼ 3%) of unobscured, scattered light, and 3) the dust emission of the surrounding torus (overall similar to the Polletta et al. 2007 torus template shown in Figure 3). This model also explains the presence of both narrow and broad emission lines be-  \nuse there is a direct sight line (albeit obscured) to the Broad Line Region, and it predicts that the red optical colors extend toward the NIR, which is roughly consistent with results showing that the handful of LRDs observed in the short-wavelength MIRI bands (F560W and F770W) continue the steep SED trend (Barro et al. 2023).  \nFigure 3 indicates that the stacked SED of the SMILES LRDs agrees well with the results from previous works based on NIRCam data: the rest-frame UV up to 0.4 µ m nicely matches the emission of an unobscured, or slightly reddened, QSO, which is consistent with the low-luminosity scattered light component of the accretion disk. Remarkably, the average SED presents a bump around the MgII emission feature, indicating that the AGN makes a non-negligible contribution to the rest-frame UV spectral range. However, the scatter around the median suggests that there might be sources other than an AGN contributing to the UV, e.g., emission from stars in the host galaxy.  \nThe rest-frame optical presents a very steep slope from 0.4 to 1 µ mconsistent with the emission of a heavily obscured supermassive black hole, here indicated with the Polletta et al. (2007) circumnuclear torus template. In addition, aided by the inclusion of the JEMS mediumband photometry, the stack exhibits statistically significant peaks at the wavelengths where we would expect strong H α +[NII]+[SII] and [OIII]+H β emission. The average LRD SED might be identified with a torus up to ∼ 0 . 7 µ m; although the torus does not have optical emission lines, the lines from the central engine can masquerade as a steeper (dustier) continuum. In summary, up to the reddest bands covered by NIRCam (F444W, F460M, and F480M filters, typically probing up to H α ), the average UV-optical SED presented in this paper roughly matches an obscured QSO spectrum plus a torus.  \nHowever, the interpretation of LRD SEDs being dominated by a pure and/or complex (with a gray attenuation law) AGN model is not supported when we consider the rest-frame near-infrared (NIR) spectral region probed by MIRI, the key addition in this paper, as first presented in Williams et al. (2023a). The MIRI data (colored circles in Figure 3) and average SED are consistently lower than the predictions of the Polletta et al. (2007) torus template. A more general evaluation can be made by applying different reddening levels to a standard unreddened AGN template. In the ∼ 0.5-3 µ m region where the MIRI measurements fall, there is general agreement on the shape of the intrinsic template, as determined by averaging the behavior of large samples (e.g., Elvis et al. 1994; Richards et al. 2006; Glikman et al. 2006; Lyu et al. 2017). Lyu & Rieke (2022b) show  \nFigure 3. Stacked LRD SEDs for sources detected by MIRI (left panel) and not (right panel), normalized at rest-frame wavelength 0.4 µ m. Black points show NIRCam fluxes for individual sources, while rainbow color points show MIRI fluxes. Arrows depict 5 σ upper limits. The average SED is shown with a magenta line (10 point averages) and its rms with a magenta shaded region. The observed average SED for LRDs is compared to 5 different templates. The orange lines show an average QSO spectrum (see text for details), and the same template extincted by A(V) = 2 mag using a Calzetti et al. (2000) attenuation law. The red lines stand for the torus template in Polletta et al. (2007), normalized to the same wavelength as the observations (continuous line) as well as normalized at the 2 µ m average SED level (dashed line) in the case of the MIRI detections plot. The blue line shows the model for an intense starburst presented in Siebenmorgen & Krugel (2007), more specifically, the sub-LIRG model with a size of 350 pc, 90% of the total luminosity coming from OB stars, optical attenuation of A(V) = 36 mag, and 10 3 cm 3 density of dust in hot spots (gas clouds surrounding and directly heated by OB stars).  \n<!-- image -->  \nthat this is an appropriate average behavior including the variations such as warm dust deficient and hot dust deficient cases (Lyu et al. 2017). It is therefore the appropriate foundation for fitting the LRD photometry, but with the large attenuations required to match the SED slope. To be specific, we adopt a SED from Glikman et al. (2006), extended to the mid-infrared (MIR) with models from Siebenmorgen et al. (2015). Figure 3 shows that there is no suitable solution; if reddening is selected to match the behavior between 0.5 and 1 µ m, the SED falls substantially higher than the MIRI photometry near 2 µ m (rest). Reducing the reddening to solve this problem yields a SED too low near 1 µ m. Overall, the MIRI photometry indicates a change to a significantly bluer slope at ∼ 1 µ m and even stabilizes at a roughly flat value at 1-2 µ m, where stellar emission is expected to peak. This is a fundamentally different shape from AGN SEDs. There are also hints of another steepening at > 3 µ m, but this behavior is traced only by a small subset of LRDs detected beyond 1.6 µ m restframe (see discussion in § 5.2); deeper MIRI data at the longest wavelengths would be necessary to probe this region.  \nThe comparison of photometry and templates in Figure 3 strongly suggests that the rest-NIR emission of the LRDs (around 1-2 µ m) is very different from the hot-dust (T ∼ 1500 K) dominated emission of the typical QSO templates used in recent LRD papers dominated by the fit of NIRCam data alone (e.g., Barro et al. 2023, Labb'e et al. 2023b, or Greene et al. 2023). Reproducing the observed MIRI colors would require a more flexible modeling of the torus emission that can vary its IR luminosity (i.e., the relative flux level with respect to other components such as stars or the accretion disk) and the wavelength at which the emission peaks. This can be accomplished using either radiative transfer models (e.g., Nenkova et al. 2008; Stalevski et al. 2012, 2016; Siebenmorgen et al. 2015), modified blackbody templates (Kim et al. 2015, Hern'an-Caballero et al. 2016 or Killi et al. 2023), or empirical templates that can account for separate contributions of the torus dust emission at different temperatures (e.g., warm and hot, as in Lyu et al. 2017). However, such models require significant departures from typical AGN behavior in the red and near infrared. Alternatively, the result could indicate that the SED at 0.5-2.0 µ m is not AGN-dominated but rather stellar-dominated, as recently discussed in  \nWilliams et al. (2023a). We demonstrate this point by plotting in Figure 3 the torus template normalized to the average SED level of LRDs at 2 µ m (dashed red line). If the emission at rest-frame wavelengths longer than 2 µ m is dominated by a torus, then the flux at around 1 µ mmust be dominated by a different component, such as the accretion disk or, more likely, stars, whose emission peaks exactly in that spectral zone (except for very young ages).  \nThe right panel of Figure 3 shows the average SED of LRDs that are not detected by MIRI. The UV range is bluer than in the MIRI-detected case, more in line with the slope of the average QSO spectrum. The change in slope at 0.4 µ m is also clear (this being one of the main selection criteria of our sample). Strong emission lines might also be present; indeed we observe a similar SED rise in the [OIII] region, and two NIRCam points may be revealing strong H α emission. However, the small number of sources does not support detecting emission lines in the average SED as well as in the case for MIRI detections. For these sources, the MIRI upper limits also imply that the SED flattens at ∼ 1 µ m, i.e., stellar emission could dominate this spectral range. However, the SED possibly has a steepening compatible with the presence of strong hot dust emission at wavelengths around 2 µ m and redward; the currently available MIRI data cannot constrain this spectral region appropriately.", '3.2. MIRI colors': 'Figure 4 illustrates further what the MIRI observations reveal about the nature of LRDs, probing further into the redshift effects. Two of the top panels show the expected NIRCam-MIRI colors for different types of templates as a function of redshift. The nature of the LRD emission in the spectral range probed by the F444W-F560W or F444W-F770W colors cannot be disentangled: both AGN and stellar dominated emission would present very similar color differences of about 0.5 mag, and the influence of emission lines in one and/or the other filter at specific redshifts complicates the problem. Only when observing at longer wavelengths with MIRI, 10-15 µ m, would the data start to differentiate between the models. For our detections at these wavelengths, the MIRI data reveal flat SEDs, F770W-F1000W, and F770W-F1280W values, bluer than 0.5 mag for 60% of the sample and consistent with stellar dominated emission presenting a 1.6 µ m bump. We note however that stellar emission alone leads to flat or even negative MIRI colors, whereas the median colors of the LRDs, including the upper limits, show small, but increasingly redder colors with wavelength which suggest that there is emission from another  \nsource, such as nebular continuum or dust, but not as dominant as in the QSO templates (see also discussion in § 5.2).  \nIn summary, as found in the preceding section, this more detailed analysis shows how SEDs constructed up to the reddest NIRCam bands are compatible with a dominant AGN with composite emission: blue in the UV, with emission lines coming from the broad-line region, and starting to rise as expected for dust torus emission up to 5 µ m. MIRI data at 5.6 µ mdeviates very little from this scenario, but at 7.7 µ m starts to tell a different story. Even longer wavelength data clearly point to a restricted contribution of a dust torus at rest-frame wavelengths around 1-2 µ m, where the stellar emission is expected to peak.  \nBoth in Figures 3 and 4, we show one of the radiative transfer models of dust-rich compact nuclear starbursts and luminous infrared galaxies (LIRGs) presented in Siebenmorgen & Krugel (2007). In particular, we plot (in blue) the template for an extreme starburst with a 350 pc star-forming region with a sub-LIRG-like luminosity (10 10 . 1 L ⊙ ; but LIRG-like luminosity density), with 90% of its luminosity arising from OB stars embedded in very thick [A(V) = 36 mag] and dense (10 3 cm -3 hydrogen density) birth clouds. Remarkably, this stellar-only model, mainly based on the presence of highly embedded OB stars within a more general stellar population, nicely reproduces the main characteristics of LRDs across the whole spectrum: the change in spectral slope at rest-frame ∼ 0 . 4 µ m and the flattening of the SED probed by MIRI. It, however, presents a steeper slope in the FUV, but we remark the stellar emission is constant across all these models and does not take into account nebular emission. The redshift dependence of the color of this template, shown in Figure 4, is also completely consistent with the measurements.  \nWith the general trends discussed with Figures 3 and 4 in mind, in the following section we investigate several different new scenarios to model the NIRCam+MIRI SEDs of each LRD in our sample.', '4. SED MODELING CODES FOR LRDS': 'In this section, we present the detailed SED fitting methods that will be applied to the analysis of individual objects in our sample to infer their physical properties. The results will be presented in the next section, concentrating on those galaxies that have spectroscopic redshifts and are detected in several MIRI bands. The SEDs of all our LRDs were fitted with 4 codes, each one with different assumptions and approaches: synthesizer-AGN , prospector-SF , prospectorAGN , and prospector-AGN +.  \n<!-- image -->  \n<!-- image -->  \nFigure 4. NIRCam and MIRI colors vs. redshift for all LRDs detected at least in F770W. The lines illustrate the color-redshift tracks for the same templates shown in Figure 3. Galaxies from different samples are marked with different colors. The first two panels in the top row show the F444W-F560W(F770W) colors which are available for the majority of the sources. The last panel at the top and the bottom panels show the MIRI-MIRI colors relative to F770W. The median and scatter of the colors are indicated in the top left corner.  \n<!-- image -->', '4.1. The synthesizer-AGN code': "For the synthesizer-AGN code (P'erez-Gonz'alez et al. 2003, 2008) we assumed a composite stellar population with both a young and a more evolved starburst. The stellar emission is described by the Bruzual & Charlot (2003) models, assuming a Chabrier (2003) IMF with stellar mass limits between 0.1 and 100 M ⊙ . Both starforming events are characterized by a Star Formation History (SFH) described by delayed exponential function with timescales between 1 Myr and 1 Gyr, and with ages from 1 Myr up to the age of the Universe at the redshift of the source. Nebular emission is also considered (see P'erez-Gonz'alez et al. 2003 for details). The attenuations of the stellar and nebular emission are described by the Calzetti et al. (2000) law, with A(V) values ranging from 0 to 4 mag for each population, considered to have completely independent attenuations. The dust emission is also modeled (beyond what the AGN might contribute) with the radiative transfer templates of nuclear starbursts presented in Siebenmorgen & Krugel (2007). We use the templates for 350 pc diameter starforming regions, with different SED shapes parametrized by the total IR luminosity (from sub-LIRG to Hyper-  \nLIRG ranges), visual extinction of the starburst (from 2 to 144 mag), ratio to the total luminosity (from 4090%) of the emission arising from OB stars embedded in dense molecular clouds of a variety of hydrogen number densities (from 10 2 to 10 4 cm 3 , assuming a gas-to-dust ratio of 150).  \nApart from the stellar population, we include an AGN component described by a QSO template constructed with the average spectrum of QSOs presented in Vanden Berk et al. (2001) and Glikman et al. (2006), extended to the far-IR with the dust torus template in Siebenmorgen et al. (2015). This model was shown in Figure 3. The attenuation of the AGN model is also assumed to follow the Calzetti et al. (2000) law, and we impose an energy balance criterion so the energy absorbed by dust is reradiated in the mid- and far-IR, scaling up the dust torus emission.  \nThere are 9 free parameters to describe the stellar population. Apart from the 4 parameters for each stellar population, i.e., age, timescale, attenuation and metallicity, we fit the stellar mass ratio between the youngest stars and the total mass, the burst strength. We add to the 9 stellar parameters, 2 more that describe the scal-  \ning of the AGN template with respect to the stars and its attenuation.  \nWe remark that synthesizer-AGN assumes independent extinctions for the 2 stellar populations in the host (as well as for the AGN). This is different from what is done by the next codes, where the attenuation for the stellar emission is governed by common parameters (visual attenuation and attenuation law slope), although internally young stars are treated differently from old stars.", '4.2. The prospector-SF code': 'For prospector-SF (Leja et al. 2017, Johnson et al. 2021), we use the MIST stellar evolutionary tracks and isochrones (Choi et al. 2016), a Chabrier (2003) IMF, a range in stellar metallicity between -1.0 and 0.19, and gas phase metallicity log(Z / Z ⊙ ) between -2.0 and 0.5. For the SFH, we use a non-parametric model, following the flexible SFH prescription (Leja et al. 2019a) with 6 time bins and the bursty-continuity prior (Tacchella et al. 2022). The ionization parameter for the nebular emission ranges from log U = -4 to -1. The nebular line and continuum emission are generated using CLOUDY (Ferland et al. 1998). We base the attenuation law on a dust model that combines: (1) the Charlot & Fall (2000) two-component approach, birth-cloud vs. diffuse dust screens; and (2) the Kriek & Conroy (2013a) method that parametrizes the diffuse component as a combination of a Calzetti attenuation plus a Lorentzian Drude to model the strength of the UV bump. Both components of the diffuse dust are then modulated by a power-law factor n that varies the slope of the attenuation between -1 to 0.4 relative to Calzetti ( n = 0).  \nThis model includes 5 free parameters that control the ratio of SFR in six adjacent time bins; the first two bins are spaced at 0 -5 Myr and 5 -10 Myr of lookback time, and the remaining four bins are log-spaced to a maximum age of 100 Myr. In addition, it fits (1) for the ratio of the nebular to diffuse attenuation, which ranges between 0 and 2, but follows a clipped normal prior centered on 1, and also (2) for the dust index n .', '4.3. The prospector-AGN code': "A third type of fit, prospector-AGN hereafter, makes use of a hybrid galaxy+AGN model similar to the one described in Barro et al. (2023). Briefly, we use a combination of a stellar emission component that dominates in the rest-frame UV and an AGN one that dominates in the optical to IR wavelengths. By construction, this is the only model where the AGN makes up the bulk of the luminosity and, as such, it serves as a feasibility test for AGN-dominated SEDs (forced to be  \ndominant, not free as in the case of the prospectorAGN + code presented below). The stellar emission is modeled with Prospector using a parametric delayedτ SFH with no attenuation. The AGN has two distinct components that account for the emission of the accretion disk and the torus, respectively. The accretion disk is modeled after the empirical QSO templates mentioned above from 0.1 µ m up to 0.6 µ m rest-frame, followed by a power-law (f ν = ν α , see for example Hern'an-Caballero et al. 2016 or more recently in Bosman et al. 2023 and Killi et al. 2023) with variable spectral index ranging between α =0 to 0.5, with default value of α = 1 / 3, and attenuated with a Calzetti et al. (2000) law. The torus emission is modeled using the clumpy torus models from Nenkova et al. (2008) included in prospector (Leja et al. 2018). These have two free parameters: the optical depth, which ranges between τ V = 1 and 150, and the total IR luminosity. Since these parameters are independent of each other, the models allow for AGN-dominated SEDs where the emission from the torus dominates at longer wavelengths ( λ ≳ 2 µ m), as opposed to the typical QSO templates where the torus emission dominates at λ ∼ 1 µ m or blueward.", '4.4. The prospector-AGN+ code': 'Last but not least, prospector-AGN + is a modified version of the original prospector code with the Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) for the stellar component. We assumed a Kroupa (2001) initial mass function and delayed-tau star-formation history. The stellar nebular line and continuum emission are also turned-on, as preconfigured in FSPS (Byler et al. 2017). We adopted the Calzetti attenuation curve with a flexible slope as introduced in Kriek & Conroy (2013b). For the galaxy dust emission, since the objects of this work have z ≳ 5 . 0, we adopted the empirical IR SED model of Haro 11, a lowmetallicity star-bursting dwarf galaxy that is believed to share typical features of first-generation galaxies in the early Universe (Lyu et al. 2016; De Rossi et al. 2018).  \nThis code uses a set of semi-empirical AGN SED models, which have been optimized for AGN identification and characterization (Lyu et al. 2022, 2023). For this work, we adopted a similar model configuration as in Lyu et al. (2023) for the SMILES+JADES AGN identification: the AGN component includes both the AGNpowered continuum from the UV to the far-IR and the narrow and broad emission lines from the UV to the NIR derived based on empirical observations. The continuum shape and line strengths of the AGN SED can be modified by a hybrid extinction configuration featuring the SMC-like curve for the commonly seen UV-optical  \nextinction in Type-1 AGNs and an empirical attenuation law for the IR obscuration. Further details can be found in Lyu et al. (2023) and references therein.  \nIn total, this code uses 7 free parameters to describe the stellar component and 2 for the AGN.', '4.5. Summary': 'These four codes let us compare models of the LRDs with a wide variety of input assumptions. This is important to let us deduce which properties are likely to be intrinsic rather than due to modeling issues. synthesizer-AGN , prospectorSF , and prospector-AGN+ use three different prescriptions for the stellar populations, with three varieties for the interstellar extinction law. synthesizerAGN and prospector-AGN+ treat the stellar-heated dust emission in two different ways. synthesizerAGN , prospector-AGN , and prospector-AGN+ use three different models for the AGN emission. This latter variety is particularly important given the diversity of possibilities for AGN SEDs and the importance of capturing this variety in modeling LRDs. The first case, synthesizer-AGN , uses an empirical fixed model SED. prospector-AGN + fits empirical templates obscured by an attenuation law derived specifically for AGNs (see e.g. Lyu et al. 2017), while prospector-AGN uses a declining power-law with α = 1 / 3 in place of the empirical templates.', '5.1. Galaxy by galaxy analysis for the Golden Five LRDs': 'In this section, we present our analysis of the SEDs of the LRDs detected in several MIRI bands. We concentrate our discussion of properties on the five galaxies detected in F1280W and redder bands, which we call the Golden Five sample. We compare the global results of the SED fitting provided by each one of our 4 different codes in different subsections, and comparatively discuss how they fit the UV, optical and near-infrared spectral regions. At the end of the section, we describe the properties for the galaxies detected up to F1000W and we briefly discuss the rest of the sample. A more detailed discussion of all the galaxies in our our sample is presented in Appendix C. The main physical properties of all the LRDs are given in Table 2.', '5.1.1. Results for the Golden Five galaxies with synthesizer-AGN': "With synthesizer-AGN , the rest-frame UV spectral range is fitted with a young stellar population for all the galaxies in the Golden Five sample except for  \nJADES -57356, which is fitted with the QSO component. For this galaxy, the fact that no emission lines are observed in the optical ([OIII] or H α ), but the UV spectral range is flat, favors the contribution of a QSO to the UV emission, outshined by the stars in the optical. Very young stars would contribute strong emission lines in the optical (as in JADES -204851), and thus are disfavored by this modeling. The different UV (relatively flat) slopes observed for the Golden Five sample are reproduced with ∼ 1 Myr old (mass-weighted age) stellar population and A(V) = 0 . 6 mag dust attenuation, which is also able to reproduce the emission lines observed for four of the Golden Five galaxies.  \nThe rest-frame optical range of the Golden Five sample is dominated by slightly more evolved stars in all cases, typically of 10-100 Myr mass-weighted age. They account for most of the stellar mass (typically, more then 90%) but do not contribute much, in relative terms with respect to the younger stars, to the optical emission lines. This more evolved stellar populations present high attenuations, typically A(V) = 2 -3 mag.  \nIf we consider the ratio of the attenuation of the farUV with respect the optical emission arising from all the stellar populations included in the synthesizerAGN modeling, κ FUV , we typically find flatter slopes, ∼ 1 . 5 -2 . 0, than what is expected for the Calzetti et al. (2000) attenuation law ( κ FUV = 2 . 6). This effect is, in fact, common to all modeling techniques, they all typically obtain gray attenuation laws.  \nFinally, the rest-frame near/mid-infrared emission of the Golden Five LRDs is fitted with dust combined with star formation in three Golden Five galaxies, and dust in an AGN torus in the other two. The heating source is indicated by the possible detection of a most probably star-formation related 3 µ m PAH feature in JADES -57356, and by the differences in slopes of the SEDs for the other sources, steeper in the case of AGNdominated fits (JADES -211388 and JADES -79803).  \nThe synthesizer-AGN code presents the smallest χ 2 values for the fits of the SEDs of all the Golden Five galaxies except one.  \n5.1.2. Results for the Golden Five galaxies with prospector-SF  \nBy construction, prospector-SF only considers stars and dust heated by stars to fit the SEDs. In general, the χ 2 values obtained by this code and the following one are larger than those obtained with synthesizer-AGN (except for JADES -219000) and smaller than prospector-AGN +.  \nThe prospector-SF results imply a ∼ 30 Myr stellar population (mass-weighted age), with large attenuations A(V) = 3 mag and a flat attenuation law. These  \nFigure 5. SED fitting results for source JADES -57356, the LRD in our sample detected up to F2100W. The four upper panels show: (1) the fits for synthesizer-AGN + on the top left panel, including the individual components of the model, i.e., young, old, and all stars in blue, cyan, and gray, (un)obscured AGN in orange, and regular dust emission (in principle, linked to star formation) in red; (2) results for Prospector-AGN + on the top right panel, showing the three components, stars in gray, AGN in orange, and star-formation heated dust in red; (3) Prospector-SF on the bottom left, showing stars in gray and star-formation heated dust in red; and (4) Prospector-AGN on the bottom right, including young stars in blue, total AGN emission in orange, with the torus emission shown with a dashed line. Number of bands fitted and direct (i.e., not reduced) χ 2 values are provided, as well a stellar masses, stellar mass-weighted ages, V -band stellar attenuation, ratio between the FUV and optical stellar attenuation, and fraction of bolometric luminosity coming from the AGN. The fits include NIRCam bands, shown in black, and MIRI fluxes, in color; upper limits, depicted with triangles, are also used. Below the SEDs, we give 10 '' × 10 '' postage stamps in NIRCam (upper row), MIRI (middle row), and MIRI convolved with a 5-pixels wide tophat filter, the LRD being marked with a 0.3 '' radius circle, and the S/N provided (when it is above 5).  \n<!-- image -->  \nFigure 6. Same as Figure 5, but for source JADES -204851, a source detected up to F1800W.  \n<!-- image -->  \nFigure 7. Same as Figure 5, but for source JADES -219000, a source detected up to F1500W.  \n<!-- image -->  \nFigure 8. Same as Figure 5, but for source JADES -211388, a source detected up to F1280W.  \n<!-- image -->  \nFigure 9. Same as Figure 5, but for source JADES -79803, a source detected up to F1280W.  \n<!-- image -->  \nproperties translate to larger stellar masses compared to other codes (a ∼ 0.5 dex difference with respect to synthesizer-AGN , for example). We remark that the dust emission models are significantly colder than what is obtained with synthesizer-AGN (see Appendix B).", '5.1.3. Results for the Golden Five galaxies with prospector-AGN': 'The main distinct characteristic of the results achieved with prospector-AGN for the Golden Five galaxies is the smaller stellar mass. Given that, by construction, this code only fits the rest-frame UV spectral region with stars, and the optical and infrared with AGN templates (including emission from the accretion disk and the torus), the stellar masses are up to a factor of 100 smaller than those estimated with the other codes. In general, better fits in terms of χ 2 values are obtained with prospector-AGN compared to prospector-SF , except for JADES -57356, the galaxy with possible PAH emission. But better fits are obtained with the mixed models in synthesizer-AGN (favoring stellar emission in the UV and optical, a variety of results in the nearinfrared).', '5.1.4. Results for the Golden Five galaxies with prospector-AGN+': 'The prospector-AGN + code fits to the SEDs of the Golden Five galaxies are very similar to those outlined for synthesizer-AGN . The rest-frame UV and optical spectral ranges are dominated by stars, except at 0.8-1.0 µ m for two galaxies (JADES -211388 and JADES -79803) whose infrared emission was fitted with an AGN torus model. This was also the best solution for synthesizer-AGN . The only significant difference with respect to the results achieved with synthesizer-AGN is found for the infrared emission of JADES -219000, whose slope if we take the MIRI upper limits at face value is too steep for a star formation model, so prospector-AGN + (which used upper limits as regular points) prefers a torus template (while the synthesizer-AGN star-formation heated dust model lies below the upper limits). The typical mass-weighted ages for the fits with the prospector-AGN + code are around 150 Myr, with extinctions around A(V) = 1 . 5 mag with a flatter attenuation law compared to the Calzetti et al. (2000) recipe.  \nA detailed analysis of the SED fits for each galaxy in the Golden Five sample is presented in Appendix C.  \n5.2. Implications for the nature of the NIR emission of the LRDs  \nOverall, all four models fit the characteristic bimodal SED of the Golden Five LRDs relatively well. Qual-  \nitatively, the UV-optical spectral range is dominated either by 1) stars: two young populations with very different attenuations in synthesizer-AGN , or a single population with an extremely gray attenuation in prospector-SF and prospector-AGN +; or 2) emission from an obscured accretion disk combined with a young stellar population that contributes only in the UV. In the next section, we discuss the implied stellar masses and other stellar properties for the Golden Five and all other LRDs that have very similar UV-tooptical SEDs (see, e.g., Figure 11) probed primarily by the NIRCam bands. However, the Golden Five, and other MIRI-detected subsamples, allow us to probe further into the rest-frame NIR of the LRDs to characterize the origin of the emission in that spectral range.  \nOur results confirm the flattening of the LRD SEDs in the rest near infrared. Previous LRD papers fitting NIRCam-only SEDs with empirical AGN templates (e.g., Kocevski et al. 2023; Barro et al. 2023; Greene et al. 2023) implied that the steep rest-optical slope would continue into the NIR. This trend seemed to be confirmed by the handful of LRDs with MIRI detections up to F770W and F1000W (Akins et al. 2023; Barro et al. 2023). However, as shown in section 3.1 and also Williams et al. (2023a), the MIRI data at longer wavelengths indicate that LRDs have a flattening in the SED between 1 and 2 µ m (rest). Interestingly, the SEDs of four out of the Golden Five LRDs at redshifts z < 7, which have direct detections of the NIR continuum at λ ≳ 2 µ m, appear to show upturns with different slopes at longer wavelengths. This suggests that, while it does not dominate, the amount of dust emission from star formation or an AGN can vary substantially from object to object.  \nOur four codes fit the NIR spectral range with a combination of the dominant source of UV-optical emission (i.e., stars or an accretion disk) plus a variable contribution from dust emission, either from star-formation ( synthesizer-AGN , and prospector-SF ) or from an AGN torus ( prospector-AGN , and prospectorAGN + depending on the source). The left panel of Figure 10 illustrates this variation showing the rest-frame color-color diagram for the LRDs and some templates and models. The 0.4-to-1 µ m color probes the optical to NIR slope. This color is similar to the F277WF444W used in the sample selection but it is not affected by emission lines and, thus, is a better proxy for the amount of dust attenuation. The color is also similar to the rest-frame V -J which is a known tracer of large dust attenuation in the UVJ diagram ( V -J > 1 . 5 mag for very dusty galaxies, e.g., Brammer et al. 2011, Wuyts et al. 2011).  \nFigure 10 illustrates how the modeling codes populate the color-color diagram between the flat (color ∼ 0 mag) stellar-only sequence, and the hot-dust dominated sequence of the QSO1 template with increasing contributions from dust emission (i.e., larger IRluminosities) from star-formation or an AGN, relative to the stellar or accretion disk continua. The color-color  \n<!-- image -->  \nFigure 10. Rest-frame optical and NIR colors for the LRDs (left, with measured and upper limit fluxes shown with circles and arrows, respectively) and different templates and best-fit models (right). The 0.4-to-1 µ m color traces the optical slope and thus it is a good proxy for the dust attenuation. The 1-to-3 µ m color tracks the amount of dust emission relative to the stellar or accretion disk emission which dominates the UV-optical SED, probed by the 0.4-to-1.0 µ m color. The LRDs exhibit values in between a stellar-only sequence with 1-to-3 µ ∼ 0 (or up to 0.5 mag with increasing nebular continuum) indicated by the solid black and magenta dashed lines on the right, and the torus-dominated sequence outlined by the colors of the Polletta et al. (2007) QSO1 template with increasing A(V), indicated with a solid grey line on the top-left, converging to the color of a Polletta et al. (2007) Torus template (black square). The dashed and dashed-dotted grey lines show a similar sequence for the HDD template of Lyu et al. (2017) and the accretion disk model with slope α = 1 / 3 used in prospector-AGN . The green and blue dashed lines and markers illustrate the sequence toward redder 1-to-3 µ m colors with increasing dust emission from star-formation (log( L IR /L ⊙ ) = 12 -13) or the Nenkova et al. (2008) clumpy torus (log( L IR /erg s -1 )=44.0-45.5) relative to the stellar or accretion disk continuum (solid black and dashed grey lines on the right).  \n<!-- image -->  \nThe behavior can be discussed in terms of the colorcolor behavior in Figure 10. The black horizontal line in the figure shows the flat 1-to-3 µ m color of a stellar population with zero contribution from dust emission, which peaks at 1.6 µ m. The magenta region indicates the redder colors up to 0.5 mag relative to the stellar continuum due to increasing amounts of nebular continuum (magenta dashed line in the left panel). The 3 grey lines indicate the color tracks with increasing attenuation (A(V)=1 to 4 mag) for the Polletta et al. (2007) QSO1 template (solid), the hot dust deficient (HDD) template of Lyu et al. (2017) used in prospectorAGN+ (dashed-dotted), and the accretion disk model with declining slope ( α = 1 / 3) and zero dust emission used in prospector-AGN . The right panel of Figure 10 illustrates some of same trends as they affect the SED templates.  \ntracks ranging from moderate to high infrared luminosities (indicated in log( L ( L ⊙ )) for the star-forming case and log( L ( ergs s -1 ) for AGNs) are computed by scaling the dust emission (f dust /f total =5% to 60% at 2 µ m) of JADES -57356. As expected, the star-forming dust templates exhibit larger luminosities than the torus at similar 1-to-3 µ m colors because their SEDs extend to longer wavelengths with a more prominent peak (see also Figure 16 in the appendix). While a single NIR color is not able to capture all the nuances of the different dust emission templates, Figure 16 shows that the torus and star-forming dust templates in Siebenmorgen & Krugel (2007) and Nenkova et al. (2008) have similar slopes for the same normalization (i.e., same f dust /f total value). Consequently, the four modeling codes can all, in principle, reproduce the NIR continuum of the LRDs.  \nThe Golden Five galaxies with direct detections beyond rest-frame ∼ 2 µ m exhibit colors that are at least 1 mag redder than a flat, stellar-only SED. In particular, two of them (JADES -57356 and JADES -79803) have very red colors, [1-to-3 µ m] ∼ 2 mag, indicative of large dust emissions and IR luminosities. JADES-219000 and JADES-211388 are not well constrained beyond 2 µ m because of their higher redshifts, but upper/lower trian-  \now the range in possible colors spanned between the prospector-SF and synthesizer-AGN best-fits, which feature different amounts of dust emission within the upper limits of the redder MIRI bands. These limits fall within the overall behavior of all the sources in the figure. The broad range overall in Figure 10 highlights the need for deep, long wavelength MIRI data to constrain precisely the amount and heating nature of dust emission in LRDs.  \nHowever, we can already see that fitting the colors with a purely AGN-dominated model requires an accretion disk model with declining slope (dashed grey line) and only a small contribution from dust torus emission (thin blue line) to successfully reproduce the bluest 1-to3 µ m ≲ 1.5 colors at the lower limit of the HDD template. This would differ significantly from lower redshift AGN, which tend to have strong emission from their circumnuclear tori. This issue is mitigated by prospectorAGN+ , which can reproduce those colors with a hybrid of AGN and stellar emission.  \nWe now discuss the LRDs detected only at shorter wavelengths. Figures 11 and 17 show the stacked SEDs for the LRDs detected up to F1000W and F770W, which reveal that their SEDs are only well constrained up to rest-frame ∼ 1.6 µ m and ∼ 1 µ m, respectively. Consequently, the 1-to-3 µ m colors of their best-fit models span a much larger range from the synthesizer-AGN , stellar-only fits with [1-to-3 µ m] ∼ 0 mag (e.g., top-left panel of Figure 11), to the much redder best-fit models of prospector-AGN+ and prospector-AGN (bottom-left and right), which sometimes fit the MIRI upper limits with pronounced upturns at λ > 2 µ m. It is worth mentioning the F770W-only sample places more restrictive constraints against very red QSO-like colors than the F1000W-only sample. This is because, as shown in Figure 4, the F1000W-only LRDs have redder F777W-F1000W colors than the upper limits of the F770W-only LRDs and thus lead to a stronger flattening of the SED in the 1-2 µ m range. Consequently, the possible upturn after 2 µ m is not nearly as red.  \nIn summary, the exact contribution from dust emission to the NIR SEDs of the LRDs is still poorly constrained by the available MIRI data at long wavelengths. Nonetheless, the constraints point to a certain diversity in the dust emission. This emission is clearly larger than stellar-only SEDs but in many cases lower than prior expectations based on QSO templates. Overall, the main conclusion is that most LRDs could harbor some (relatively large) amount of dust emission and the heating source must be intense, although not necessarily a dominant obscured AGN.', '5.3. Physical properties of the LRDs detected in the bluest MIRI bands': 'All other LRDs in our sample detected up to F1000W, apart from the Golden Five, are shown jointly in the SED plots provided in Figure 11. In the same way, the SED fits for all F770W detections not included in any previous SED plot, as well as all sources not detected by MIRI are shown and discussed in Appendix C. Physical properties derived from each code for individual galaxies as well as statistical properties obtained for the whole sample and subsamples are given in Tables 2 and 3.  \nOverall, the properties of those samples with fewer MIRI points in their SEDs are similar to the Golden Five galaxies; the clearest difference is that some blue and/or very flat SED sources start to enter the selection. These are selected due to a very strong F444W emission probably linked to a high-EW emission line, but the slope of the SED is not very different in the rest of the SW and LW filters, with a quite flat slope (or even blue, as in the case of JADES -187025).  \nEven though the dust emission spectral region is not fully probed for the 10 µ m sample, and to an even lesser extent the 7.7 µ m and non-MIRI samples, the upper limits imposed by the MIRI data at longer wavelengths, more specifically, at 12.8 and 15 µ m, indicate a very similar behavior of the spectral range around 2-3 µ m (rest) compared to what we showed for the Golden Five galaxies. Indeed, the SED flattens, indicating that possible dust emission powered by stellar or AGN heating is not dominant and could only start adding significant flux redward of 3 µ m. The properties we infer from the UVto-NIR SEDs are also similar for the F1000W, F770W, and no-MIRI subsamples compared to those obtained for the Golden Five galaxies. There are, however, some observational trends, which translate to differences in physical properties.  \nFirst, the Golden Five galaxies are brighter than the rest of the sources in the full sample (cf. Figure 1). The median and quartiles for F444W are 25 . 4 26 . 1 25 . 1 mag, compared to 26 . 1 26 . 3 25 . 8 mag for the F1000W sample, 25 . 9 26 . 1 25 . 7 mag for the F770W sample, and 27 . 6 27 . 8 26 . 8 for the sources with no MIRI detection. Concerning colors, the F1000W sample is redder than the Golden Five, and the overall sample. As shown in Figures 1 and 4, the F1000W sample is among the reddest in both the NIRCam and MIRI colors: F277W-F444W=2.2 mag, F444W-F770W=1 mag and F770W-F1000W=0.9 mag, versus for the overall sample they are 1.4 mag, 0.7 mag and 0.2 mag respectively. Comparatively, the largest difference is in F770W-F1000W color difference, where the F1000W sample has similar colors to 2 galaxies among the Golden Five, JADES-79803 and JADES-  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 11. SED fitting results for the 7 galaxies detected at F1000W and not beyond (i.e., the plot does not include any of the Golden Five galaxies). SEDs are normalized to 0.7 µ m. In gray, we show the fits to each individual galaxy, and we provide an average in brown.  \n<!-- image -->  \n211388, whose MIR spectral range is fitted with dust tori. We note that for JADES-211388, at z sp = 8 . 3846, F1000W lies on top of the He I+Paγ line, which partially explains the red color.  \nLooking at the rest-frame colors and stacked SEDs in Figures 10 and 11, we find that the F1000W sample has also redder 0.4-to-1.0 µ m colors than the Golden Five sources, suggesting that they are dustier (see next paragraphs). Interestingly, using a longer baseline color 0.25-to-1.0 µ m (similar to NUV-J), which probes into the relatively flat UV SED of the LRDs, we find even larger differences between the Golden Five galaxies and the F1000W samples, 0.25-to-1.0 µ m = 2.1 mag for the Golden Five vs. 3.5 mag. These colors and the multicomponent SED modeling discussed in the previous section indicate that the colors relative to rest 1 µ m are partially driven by differences in the relative luminosity  \nof the component dominating the rest-UV (young, unobscured stellar population) and the component dominating the rest-optical (older stellar population or obscured accretion disk). As the flat UV component scales up in brightness, it leads to bluer 0.25-to-1.0 µ m and 0.4to-1.0 µ m colors, but perhaps not because of a change in the intrinsic properties of the component dominating the optical range. It could be that a brighter UV reveals a larger fraction of the starburst emission percolating through the compact dust cloud (possibly linked to a higher burst strength or younger age, apart from dust-star relative geometry), or perhaps it shows a more massive stellar host for the obscured AGN. This interpretation also helps to explain the larger scatter in the UV region of the stacked SEDs relative to the optical region. That is, while all LRDs have distinctive blueUV and red-optical SEDs, there is a larger diversity in  \nthe UV emission for a similar optical-to-NIR slope that might reflect variations in the relative luminosity of two different components.  \nThe overall colors of the LRDs range from 0.25-to1.0 µ m = 2 to 4.5 mag. Taking 0.25-to-1 µ m = 3 mag as an intermediate value, we find that all of the Golden Five galaxies exhibit bluer colors, versus only 30% (2/7) of the galaxies in the F1000W and F770W samples. This emphasizes again that the Golden Five galaxies are intrinsically bluer than the other samples in all colors. Interestingly, we also find a trend toward stronger emission lines (larger EWs) with bluer 0.25-to-1.0 µ mcolors. For example, this trend is seen among the bluest Golden Five sources (JADES -79803 and JADES -204851) and in the handful of galaxies with bluer colors in the F1000W sample (JADES -210600, JADES -214552 and JADES -217926) and F770W sample (JADES -187025 and JADES -197348).  \nBased on the previous observational differences, our modeling of the SEDs provides some trends also in physical properties. The Golden Five galaxies lie at smaller redshifts than the rest of sources, median values z = 5 . 5, compared to z = 6 . 3, z = 5 . 9, and z = 7 . 4 for the F1000W, F770W, and non-MIRI subsamples (check Table 3 for more statistical information). So part of the reason for the detections of the Golden Five in many MIRI bands can be linked to redshift.  \nThe stellar masses of the Golden Five galaxies are 0.20.4 dex larger than those of the other subsamples, with values around 10 9 . 6 M ⊙ for the former. Interestingly, the stellar population attenuation is lower for the Golden Five galaxies, A(V) = 3 . 0 mag, compared to F1000W sources (A(V) = 3 . 2 mag), F770W galaxies (A(V) = 3 . 3 mag), and non-MIRI sources (A(V) = 3 . 9 mag).  \nWe conclude that sources detected only in the bluest MIRI filters are not just fainter (less massive) versions of the Golden Five galaxies, lying at higher redshifts. Indeed, there are also other differences in physical properties (deriving from differences in SEDs) which point to larger attenuations. If the fainter MIRI sources present larger attenuations, but still are not detected by MIRI at the longest wavelengths, the interpretation would be that the host dust emission is not enhanced compared to the Golden Five galaxies, which would lead to a dominant role of dust-enshrouded star formation rather than obscured nuclear activity for a significant fraction of these sub-populations.', '5.4. Statistical stellar and AGN properties of LRDs': 'In this section, we discuss the general properties of our sample of LRDs. Figure 12 shows the distributions of stellar mass, bolometric luminosity, mass-weighted age,  \nand bolometric stellar light attenuation. All statistical information is summarized in Table 3.  \nThe typical stellar mass of LRDs is log M ⋆ / M ⊙ = 9 . 4 9 . 7 9 . 1 (median and quartiles) according to synthesizer-AGN . Considering the full redshift range of our sample of LRDs, 5 ≲ z ≲ 9, and the number density of galaxies in the stellar mass range 9 . 0 < log M / M ⊙ < 10 detected by CEERS based on the v0.51 catalogs (Finkelstein et al. 2023), we calculate that the LRDs account for 14 ± 3% of the full population of galaxies (subject to uncertainties due to cosmic variance). This translates to a comoving density of LRDs of 10 -4 . 0 ± 0 . 1 Mpc -3 , which is quite constant across the 5 < z < 9 redshift range, with differences < 0 . 1 dex between 5 < z < 7 and 7 < z < 9. The estimates for both the LRDs and the other galaxies are subject to a number of potential systematic errors, but the estimate indicates that the LRDs represent a significant, but not a dominant, population over this redshift range. Given the ranges of Universe ages probed by our sample, 200 Myr for 7 < z < 9 and 400 Myr for 5 < z < 7, the ∼ 10% frequency within the global population could be interpreted in terms of a duty cycle around 20-40 Myr, which points to starburst behavior.  \nLarger masses are obtained by prospector-SF , log M ⋆ / M ⊙ = 9 . 9 10 . 2 9 . 2 , mainly because the mass-weighted ages are older, between 10 and 100 Myr, with median and quartiles being t m -w = 16 27 10 Myr. synthesizerAGN fits the SEDs with significantly younger stellar populations, typically t m -w = 3 20 2 Myr. The 2 stellar populations considered by synthesizer-AGN are typically younger than 20 Myr. In fact, the average SFH of LRDs obtained by synthesizer-AGN and shown in Figure 13 indicates that LRDs are experiencing a very intense episode of star formation extending for nearly 10 Myr and with a very compact size. The burst would be in part heavily dust-enshrouded, with some younger stars having cleared the interstellar medium and being directly observable through much smaller dust optical paths. The young ages are expected for starbursts with strong emission lines (as present in some of our galaxies), with large amounts of dust (confirmed for a significant fraction of the whole sample), and with gas also feeding a SMBH. synthesizer-SF provides a similar average SFH, extended almost at a constant level up to approximately 10 Myr, and decaying afterwards. However, the first age bin considered by synthesizer-SF encompasses the 2 bursts obtained by synthesizer-AGN , the former adding more mass in ages around 10 Myr (with a large scatter, shown by the shaded region in Figure 13).  \nComing back to stellar content, even smaller masses compared to synthesizer-AGN are obtained by  \nFigure 12. Statistical stellar properties of LRDs, according to the 4 SED-fitting codes described in Section 4. From top to bottom, left to right, we show stellar masses, bolometric luminosities (obtained by integrating the stellar emission correcting for the effects of dust attenuation), mass-weighted ages, and bolometric stellar luminosity attenuation. Medians and quartiles are shown for each distribution. For the results provided by synthesizer-AGN , we separate statistics for the young and old stellar populations (marked as you and old ) as well as the integrated values.  \n<!-- image -->  \nprospector-AGN , log M ⋆ / M ⊙ = 7 . 8 8 . 0 7 . 2 , which only considers contributions of stellar light to the UV spectral region. These masses are similar to the values obtained for the youngest population in the synthesizerAGN modeling. Given that prospector-AGN fits the optical and NIR spectral regions with an AGN, the stellar mass estimates are significantly smaller than what is needed to reproduce the emission at those wavelengths with stars.  \nWe conclude that the prospector-AGN stellar masses should be considered lower limits, since they assume little contribution of stellar light to the optical spectral range. synthesizer-AGN and prospectorAGN + typically obtain fits for which the optical emission is dominated by stars, hence the stellar mass estimates they obtain should be interpreted as more realistic or upper limits. The differences in stellar masses derived with synthesizer-AGN and prospector-SF exemplify the effect of the SFH, which is relatively important  \n(0.4-0.6 dex) for these young galaxies whose mass-tolight ratios can change significantly as the most massive stars disappear.  \nThe bolometric luminosities (including dust-absorbed energy) vary by a factor of 10 between the prospectorSF and prospector-AGN runs, while synthesizerAGN lie in between, with a typical value of log L / L ⊙ bol = 12 . 0 12 . 4 11 . 6 . Combining with the stellar masses, we conclude that LRDs present mass-to-light ratios of 1/400, typical of OB stellar associations (a B2 star would have that value, approximately), as expected based on the young ages. We remind the reader that given these very young ages and the starburst nature of LRDs, the a priori assumption of a universal IMF is quite relevant. The amount of OB stars formed, the quick and efficient formation of metals and dust, and the inferred stellar masses (or even the growth of a SMBH) are all affected by the IMF.  \nFigure 13. Average SFHs of LRDs according to the fitting codes presented in Section 4. Averages and scatter are shown as lines and shaded regions.  \n<!-- image -->  \nFinally, the bottom-right panel of Figure 12 shows the total attenuation of the stellar light in LRDs. All codes are consistent in assigning large dust content to this type of galaxy, with attenuation around 3-4 mag (i.e., 95-99% of the light being absorbed by dust).  \nThe relative importance of the AGN and stellar components in LRDs is presented in Figure 14. Here we show the fraction of the total luminosity coming from the AGN and integrated in several spectral ranges. We show bolometric, UV (integrated up to 0.4 µ m), optical (from 0.4 to 2 µ m), and IR (from 2 µ m redwards) luminosity ratios for the AGN emission (with the rest coming directly from stars or dust heated by stars). We note that we only include 3 of the 4 codes in this plot, since prospector-SF does not consider any AGN contribution (although the dust models imply intense radiation fields which could be easily identified with an AGN).  \nFor synthesizer-AGN , the bolometric luminosity of most of the sample is dominated by stars, with only ∼ 20% of sources presenting AGN luminosity fractions larger than f AGN L bol > 0 . 5. In contrast, prospectorAGN + obtains a much larger AGN contribution for most galaxies, with nearly 70% of the galaxies presenting a bolometric luminosity fraction > 50%. This is a direct consequence, however, of this code fitting the MIRI nondetections assuming the 5 σ upper limit as an actual flux. This means that the prospector-AGN + results about the AGN luminosity ratio should be regarded as upper limits. In the same sense, synthesizer-AGN could be regarded as lower limits, since the MIR is not fitted for MIRI non-detections nor for F770W-only sources, since all bands with measured fluxes are well reproduced by  \nstellar models alone (i.e., the dust emission is loosely constrained and it is not fitted).  \nMore consistent results are found for the UV luminosity fraction identified with an AGN. All codes agree that the UV spectral region of LRDs is best-reproduced by a young stellar population with varying but low dust content and relatively young ages (1-10 Myr). We remark here that prospector-AGN forces the UV to be dominated by unobscured stars, while the other 2 codes leave complete freedom in this spectral range in terms of a possible contribution from a UV-bright AGN or the amount of dust.  \nA similar behavior is observed in the optical, where both synthesizer-AGN and prospector-AGN + find that the emission is dominated by stars in most ( > 75%) of the sample. The a priori assumption of prospector-AGN is that the AGN dominates the IR spectral region and results in the optical also being dominated by the nuclear activity for most sources; consequently, as we mentioned, the estimated stellar masses are considerably smaller than the values obtained by the other codes.  \nFinally, the bottom-right panel of Figure 14 shows the fraction of the IR (dust) luminosity linked to the AGN. prospector-AGN +and synthesizer-AGN show opposite distributions, but two facts must be taken into account to interpret this behavior. First, as mentioned earlier, prospector-AGN + fits upper limits as regular flux points. Those upper limits increase with a slope similar to what can be expected for an AGN torus. Second, the dust emission fits provided by synthesizerAGN imply very intense radiation fields in dense, compact hot spots, whose properties would be indistinguishable between OB associations or an AGN (see discussion on Figure 5 and Appendix B).', '5.5. Spectroscopic properties of LRDs found in the literature': 'Apart from JADES -204851, which has broad H α emission at z = 5 . 4790 ( ∼ 2000 km s -1 ) (Matthee et al. 2023b), as mentioned in Section 5, there are other sources in our sample with relevant spectroscopic information. JADES -197348 was included in the JADES NIRSpec initial data release (Bunker et al. 2023) and identified with a broad-line AGN (Maiolino et al. 2023b). Its spectrum shows a ∼ 2500 km s -1 wide component that accounts for two thirds of the H α total flux, while [OIII] shows no such component. Our fits to this source are dominated by a QSO-like spectrum in the optical and NIR in the case of synthesizer-AGN and prospector-AGN . JADES -154428 is found to present a broad-line component with FWHM ∼ 1800 km s -1  \nFigure 14. Histograms of the fraction of integrated luminosity coming from an AGN, according to the SED fits presented in Section 4. On the top-left, we show results for the bolometric luminosity. On the top-right, for the UV luminosity, integrated up to 0.4 µ m. On the bottom-left, results for the optical luminosity are provided, with L opt defined as the integral between 0.4 and 2 µ m. The bottom-right panel shows the histograms for IR wavelengths longer than 2 µ m. Medians and quartiles for each spectral range and code are displayed at the top of the panels.  \n<!-- image -->  \n(Sun et al. 2024, in prep.); our fits include a nonnegligible contribution from a QSO-like spectrum, dominating the SED ( prospector-AGN and prospectorAGN + or accounting for ∼ 50% of the emission at specific wavelengths ( synthesizer-AGN ). No other broad H α or H β line component has been reported in the remaining 18 galaxies with available spectroscopy, i.e., 17% of the spectroscopic sample are confirmed AGN hosts. For the rest of the spectroscopic sample, the presence of an AGN cannot be ruled out, since the broad-line region could be hidden due to geometrical effects. We also note that the typical 5 σ depth of FRESCO NIRCam 3.9-5.0µ m grism spectroscopy is ∼ 5 × 10 -18 erg s -1 cm -2 for broad emission line (FWHM ∼ 1000 kms -1 ) from a point source and, therefore, faint broad line emission (possibly from an AGN not dominating the continuum) can remain undetected.', '6. SUMMARY AND CONCLUSIONS': "We characterize the nature of Little Red Dots (LRDs) in the JADES field by analyzing their spectral energy distributions including the mid-infrared fluxes provided by the SMILES program for all MIRI broad-band filters. These data probe the rest-frame near-and mid-infrared where stellar emission and/or obscured AGN emission peak. After removing brown dwarfs, which contaminate our sample at the 15% level, we arrive at a sample of 31 LRDs, the surface density being 0.9 arcmin -2 . This translates to a number density 10 -4 . 0 ± 0 . 1 Mpc -3 , accounting for 10-15% of the global population of galaxies with similar redshifts ( z ∼ 7) and stellar mass (log M ⋆ / M ⊙ = 9 . 5). Two thirds are detected in the F560W and F770W filters (all sources brighter than F444W < 26.5 mag), two fifths in F1000W, one seventh in F1500W, one thirteenth in F1800W, and one source in F2100W, down to 5 σ limits between 26.1 and 22.6 mag. The MIRI detection fraction is largely dependent on the F444W brightness, but we find an additional trend to-  \nre detections at 10 µ m or beyond with bluer F277W-F444W colors.  \nWe find that the observed MIRI colors of the LRDs, in combination with the reddest NIRCam bands, are bluer than the typical obscured QSO templates, which are dominated by the torus warm/hot dust emission at λ ≳ 1 µ m. Indeed, the rest-frame NIR spectral range exhibits a much shallower slope that is consistent with the peak of stellar emission at around ∼ 1 . 6 µ m.  \nWe modeled the rest-frame ultraviolet to mid-infrared spectral energy distributions with a battery of codes that include AGN and stellar emission templates. The various outputs allow us to identify the best fits to the distinctive short wavelength blue plus long wavelength red colors of LRDs under a range of assumptions. They also let us examine which of our conclusions are the most robust (e.g., are reflected by a number of the modeling techniques).  \nIn general, stellar-dominated models obtain a better agreement to the near-infrared at 1-2 µ mas well as in the UV at λ < 0 . 4 µ m. The AGN-dominated models where emission from an obscured accretion disk dominates the optical and NIR and the torus takes over at longer wavelengths λ ≳ 2 µ m provide a better agreement than the typical QSO templates, but still worse than the stellar models. Furthermore, this AGN-dominated model also has conceptual problems given that many of the LRDs do not present emission lines, which should be expected for the direct detection of the accretion disk (and its broad- and narrow-line regions). Consequently, we favor the interpretation that the UV-to-optical spectral range of most LRDs is dominated ( > 50% luminosity ratio) by stars.  \nIn the rest-frame near-infrared, we find that the LRDs detected in the reddest MIRI bands, beyond rest wavelengths λ ≳ 2 µ m, have color differences redder than can be expected for stars alone (even accounting for nebular emission associated with a young starburst), and consistent with some amount of emission from dust heated by star formation or an AGN. The upper limit of the MIRI colors for most LRDs also rules out that the nearinfrared emission is strongly AGN-dominated, but the loose constraints in some of them can still accommodate similar amounts of dust emission as in the LRDs detected in the rest-frame mid-infrared.  \nGiven that the rest-frame UV/optical spectral range is dominated by stars, we estimate stellar masses. The modeling of the stellar emission must consider the large attenuations implied by the red NIRCam longwavelength colors and MIRI fluxes and the need for two distinct (in terms of age) stellar populations, with differential attenuation levels or an extremely gray atten-  \nuation law, to account for the blue/flat NIRCam-SW emission as well the presence of emission lines. With this in mind, we obtain typical stellar masses for LRDs around log M ⋆ / M ⊙ = 9 . 4 9 . 7 9 . 1 , significantly smaller than what can be obtained with simple recipes for the Star Formation History (e.g., single exponential burst) and attenuation law (e.g., single Calzetti law, as noted by Barro et al. 2023). This mass estimate can be biased due to uncertainties in the Star Formation History at the 0.5 dex level (overestimated, if, for example, the stellar emission only contributes significantly to the ultraviolet emission), and can be affected by a possible (but less favored by most of our models) significant contribution of an AGN to the optical and near-infrared spectral range at the 1.5 dex level.  \nVery young stellar ages (typically around 10 Myr or younger) are supported by high equivalent width lines seen in some LRDs (with spectroscopic and imaging data), and the presence of large amounts of dust, which is a common feature of all models for the full sample, and can be expected for gas-rich, quickly enriched dense starbursts with large amounts of OB stars (affected by the IMF). However, despite the presence of large amounts of dust, LRDs are characterized by a relatively blue emission at wavelengths < 0 . 4 µ m, which can only be reproduced if this spectral range is dominated by (1) a QSO; (2) has holes in the interstellar medium surrounding star-forming regions that allow us to see unobscured very young star formation; or (3) have a gray attenuation law typically linked to significant scattering (with implications to the relative geometry of dust and stars and the dust clumpiness). Indeed, the large ∼ 3 mag scatter in the UV-to-NIR colors within the typical blue-red SED of the LRDs might indicate different burst strengths or different fractions of percolating light through the compact dust envelope. The presence of a QSO is hinted by some medium-band colors consistent with MgII emission. Nevertheless, for most of the sources, the best spectral energy distribution fit is obtained with just stars presenting a differentiated and/or gray dust attenuation law. A significant fraction of the total stellar emission of LRDs comes from OB stellar clumps mostly embedded in dense dusty regions with large optical paths, A(V) ≳ 10 mag, with the integrated stellar emission obscured at the 90-95% level, which would explain the near- and mid-infrared emission jointly with some contribution from an obscured AGN.  \nWe thank the referee for their constructive comments to our original manuscript. PGP-G acknowledges support from grant PID2022-139567NBI00 funded by Spanish Ministerio de Ciencia e Innovaci'on MCIN/AEI/10.13039/501100011033, FEDER Una manera de hacer Europa . This work was supported by NASA grants NNX13AD82G and 1255094. The work was also supported by NIRCam Development Contract NAS5-02105 from NASA Goddard Space Flight Center to the University of Arizona. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1207. DP acknowledges support by the Huo Family Foundation through a P.C. Ho PhD Studentship. BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015, and JWST Program 3215. AJB acknowledges funding from the 'FirstGalaxies' Advanced Grant from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 789056).  \nFacilities: JWST (NIRCam), JWST (MIRI), HST (ACS), HST (WFC3). All the JWST data used in this paper can be found in MAST: 10.17909/jmxm-1695 and 10.17909/8tdj-8n28.  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018), Cloudy (Ferland et al. 2013), prospector (Leja et al. 2017, 2019b; Johnson et al. 2021), synthesizer (P'erez-Gonz'alez et al. 2003, 2008).  \nTable 1. Sample of LRDs  \nNote -Table with basic information (used for selection) about the sample of galaxies in this paper: IAU-format name, ID based on JADES DR2 catalog (Eisenstein et al. 2023b), magnitude and colors used in the selection (see Figure 1), and subsample according to the reddest band counting with a MIRI detection.  \nTable 2. Physical properties of LRDs (complete version online)  \nNote -Table with physical properties of the sample of galaxies in this paper. Apart from redshift (spectroscopic values given with 4 decimals), we provide results for each of the 4 stellar population synthesis codes plus AGN described in Section 4: synthesizer-AGN (sy-AGN), prospector-AGN + (pr-AGN+), prospector-SF (pr-SF), and prospector-AGN (pr-AGN). We quote stellar masses (surviving stars), attenuation of the bolometric, far-ultraviolet (FUV, i.e., 150 nm), and visual (V, i.e., 550 nm) emission from stars (including nebular emission, and providing results for the old and young stellar populations for synthesizer-AGN in the FUV and V cases), fraction of the luminosity associated with the AGN for the bolometric, UV, optical, and IR emission (check definition in Section 5.4), and mass-weighted ages.  \nTable 3. Statistical properties of LRDs  \nNote -Table with statistical physical properties of the sample of galaxies in this paper. We provide medians and quartiles for the physical properties mentioned in Table 2 for the whole sample and the subsamples defined based on the detection in different MIRI bands.", 'A. BESPOKE PROCEDURES FOR THE REDUCTION OF MIRI DATA: THE SMILES CASE': "The MIRI data used in this paper were reduced with the Rainbow JWST pipeline developed within the European Consortium MIRI GTO Team to deal with MIRI, NIRCam, and NIRISS imaging data. The pipeline relies on the jwst official pipeline and adds some offline steps to improve the results, mainly dealing with the varying, full-of-structure background observed in MIRI data, especially at the shortest wavelengths.  \nThe Rainbow pipeline starts with a default execution of the 3 stages of the jwst official pipeline, which provides a first full mosaic, now also dealing (to some extent, but not completely) with cosmic ray showers (or snow balls for NIRCam). This mosaic is used to detect sources with sextractor in order to produce a mask for further refinements of the calibration. We use a relatively shallow detection limit since the mosaic presents intense background gradients and structure, indicating that a deep detection would be dominated by background, not real sources. Apart from the sextractor detection, we performed a 5-pixel dilation of the segmentation maps to account for the faint outskirts of extended objects, also including the emission from PSF spikes and the cruciform feature of short wavelength MIRI data (G'asp'ar et al. 2021).  \nAfter masking sources from the previous task, stage 2 data (i.e., cal.fits files) are median-filtered in rows and columns and a smooth 4-order surface is subtracted. The new calibrated data products are again mosaicked with the official pipeline. This new mosaic presents a much more well-behaved background and allows a more aggressive source detection. The new mask is used in a final step of the Rainbow pipeline which implements a super-background strategy to obtain the best results.  \nThe super-background strategy consists in homogenizing the background of a given stage 2 image using all the other images taken in the same program. For programs extending over several epochs and with enough data, PRIMER for example, we use the closest data in time. We first subtract the median of the background (after masking sources detected in the previous step) for each image. Then, for each frame to reduce, we build a stacked median background image with the rest of images. If not enough data are available for a given pixel, i.e., when very few images are available and the dithering is small compared to the size of (some) objects in the field, we replace the pixel value by a random nearby background pixel chosen from the 100 closest non-masked pixels. The stacked median super-background image is subtracted from the frame we are considering, and the filtering in rows and columns is performed, as well as the subtraction of a smooth surface.  \nThe astrometry of the new background-homogenized cal.fits files is calibrated with the tweakreg external routine provided by the CEERS collaboration (Bagley et al. 2023), using an external catalog constructed with IRAF's center task in centroid mode (which was checked to perform better than photutils), using in the SMILES case the JADES catalog as the WCS reference. Then stage 3 of the pipeline is executed, switching off the tweakreg step and setting the pixel scale to 60 milliarcsec in the case of MIRI data. The WCS of the final mosaic is checked again against the reference WCS catalog and a final background subtraction is performed with sextractor .  \nThe procedure is evaluated in Figure 15, where we compare the histograms of pixel signal for the final mosaic of the F1500W filter reduced just with the jwst pipeline (after subtracting the median value of the full image) and reduced with our super-background method. The histograms are fitted to a gaussian, showing that our procedure is able to reduce the noise by a factor of ∼ 1 . 5. This translates to 0.5 mag deeper detection limits in the final mosaic for this band compared to what can be achieved with the official pipeline alone, which roughly agrees with the expectations provided by the ETC version v3.0.  \nTable 4 .  \nNote -(1) MIRI filter. (2)-(4) Observational strategy in the SMILES survey: number of groups, number of integrations, number of dithering positions. (5) Exposure time per pixel (seconds). (6) Average background measured in the data. In parenthesis, the ETC v3.0 predictions for Dec 7, 2022 are given. Units (MJy/sr). (7) FWHM of the PSF in each filter (arcsec). (8) Depth of the data for a point-like source measure in a circular aperture with radius equal to the FWHM of the PSF, aperture corrected. In parenthesis, the ETC v3.0 predictions for Dec 7, 2022 are given. Units are AB mag. (9) Aperture correction applied to previous column results (AB mag), based on files released in pmap-.  \nFigure 15. Histogram of pixel values (transformed to µ Jy) of the F1500W SMILES mosaic, after subtracting the median background. In red, we show the results for the mosaic produced by the official pipeline. In green, the results obtained with the bespoke version of the pipeline, embedded in the Rainbow database, and implementing a super-background strategy. The 2 histograms are fitted to a gaussians, whose dispersion is translated to 5 σ depths for a point-like source, calculations based on measurements in an circular of radius equal to the FWHM of the PSF, and corrected for the limited size of the aperture using the calibration available in pmap 1138.  \n<!-- image -->  \nTable 4 provides the observation strategy, total exposure times per pixel, average background levels, and depths of the SMILES data reduced with the Rainbow JWST pipeline, compared with ETC predictions.", 'B. DUST EMISSION MODELING': "In this Appendix, we describe in detail the dust emission models used in this paper. Different recipes and emission origins are considered by the four codes described in Section 4. The main characteristics of our dust emission templates are summarized in Figure 16.  \nThe synthesizer-AGN code uses the radiative transfer models of starburst nuclei and (ultra)luminous infrared galaxies (LIRG) presented in Siebenmorgen & Krugel (2007). These models assume an intense star formation event (on top a more evolved stellar population) where a fraction of the most massive (OB) stars are embedded in compact dusty clouds (hot spots in their terminology) that dominate the mid-infrared emission. The models are parametrized in terms of the total radiated luminosity (ranging from sub-LIRG to hyper-LIRG values), the size of the star-forming region, the total amount of dust described by the total V -band attenuation, the fraction of the total luminosity linked to the OB stars in the hot spots, and the gas/dust density of the hot spot clouds.  \n<!-- image -->  \nFigure 16. Dust emission models used in the analysis of LRDs presented in this paper. Five different models from Nenkova et al. (2008) used by prospector-AGN are shown in red colors, showcasing the shift to shorter wavelengths (from ∼ 30 to ∼ 10 µ m) of the dust emission peak arising from an AGN torus with different dust optical depths, namely, τ V equal to 5,10, 20, 40, and 150. Gray lines show the dust torus and QSO templates (the latter, extincted by 1 and 2 mag following a Calzetti et al. 2000 law) in Polletta et al. (2007), the former being used by synthesizer-AGN . Thick lines show emission from dust, thin lines show the full models, which include stellar emission. Two black-body models for warm dust (500-1200 K) are depicted. The dust emission models for nuclear starbursts, also used by synthesizer-AGN (Siebenmorgen & Krugel 2007), are shown in blue and cyan colors, with representative values of the different parameters (total luminosity, luminosity arising from OB stars, hot spot hydrogen density, and total attenuation in the V -band). The star-forming model for Haro 11 used by prospector-AGN + (Lyu et al. 2016; De Rossi et al. 2018) is shown in purple, and the models from Draine & Li (2007) used by prospector-SF in green.  \n<!-- image -->  \nIn Figure 16, we show models for a 10 10 L ⊙ and 10 12 L ⊙ compact (350 pc in size) region, attenuations A(V) ∼ 10 , 20 , 100 mag, OB luminosity ratios 40 and 90%, and hydrogen number densities 10 2 and 10 4 cm -3 . The templates include a fixed stellar population, which is removed for our modeling of LRDs. Models with and without stars are shown in Figure 16.  \nWe remark that these radiative transfer models nicely recover the blue+red nature of the UV-to-NIR emission of LRDs. The change in slope of the stellar UV and optical emission is governed by the OB luminosity ratio, i.e., cyan models (OB90) change in slope at longer wavelengths, around 0.4 µ m rest-frame, similarly to LRDs. Their dust emission peaks at shorter wavelengths, implying a significant amount of warm dust. The attenuation of the FUV emission in these compact starburst models is, however, steeper than what is observed for our galaxies. All Siebenmorgen & Krugel (2007) models, which assume Milky Way type dust, present more or less prominent PAH bands. The typical mid- to far-IR flux density ratio is nearly 1000, which means that some of our galaxies, which present MIRI short wavelength fluxes around 1 µ Jy, are also constrained by the non-detections in Herschel bands (with 5 σ limits around 2-3 mJy at 100 and 160 µ m and 10 mJy at 250-500 µ m).  \nThe MIR slopes of the Siebenmorgen & Krugel (2007) models are very similar to the slopes of AGN torus emission presented in Nenkova et al. (2008), which are used in our prospector-AGN fits. Nenkova et al. (2008) parametrizes the emission from clumpy dust tori in terms of the number of clouds intercepting the visual, which depend on other parameters such as the torus thickness (outer to inner radii), the density and angular distribution of clumps and the viewing angle. The subset of templates included in FSPS and then Prospector adopts typical assumptions for the majority of these parameters (see e.g., Leja et al. 2018) leaving e only the scaling (overall luminosity) and the optical depth of an individual dust clump at 5500 ˚ A( τ V from 5 to 150) as free parameters. For the prospector-AGN fits we remove any contribution from the accretion disk in the dust tori templates and we instead model the disk separately following a combination of empirical QSO template plus a power-law f ν = ν α (equivalent to ν f ν = ν α +1 ) with variable slope α = -0 . 5 to 0.5, as discussed in Hern'an-Caballero et al. 2016 or more recently in Bosman et al. 2023, attenuated by a Calzetti law. For reference, the accretion disk model with A(V)=1.5 is shown in Figure 16 in magenta. The main difference between the Nenkova et al. (2008) models and the ones in Siebenmorgen & Krugel (2007) dominated by star  \nformation is the absence of PAH emission (especially relevant for our observations, the one at 3 µ m) and the silicate absorption present in some dust models. Overall, a dust torus (and an AGN template) is almost featureless.  \nFor reference, Figure 16 also shows the torus (also used in our synthesizer and prospector-AGN + fits) and QSO template presented in Polletta et al. (2007), the latter additionally attenuated with a Calzetti et al. (2000) law assuming A(V) = 2 mag. These templates present a similar slope as the Nenkova et al. (2008) torus models, adding the contribution from the accretion disk that dominates at wavelengths shorter than ∼ 1 µ m and a tail of colder dust peaks at around 100 µ m, as the star-forming models.  \nThe AGN SED models used in prospector-AGN + are largely based on empirical observations with the AGN emission strength calibrated against various observations. For example, the AGN hot dust emission predicted from the SED template is confirmed with NIR image decomposition of HST observations of lowz quasars (Lyu et al. 2017). The AGN mid- to far-IR SED shape is checked against MIR spectral decomposition and PAH strengths (e.g., Lyu & Rieke 2017). Compared to radiation transfer models, these empirical SED models can provide more realistic descriptions of the observations across a very wide range of AGN luminosity and redshift with fewer free parameters and less model degeneracy, which make them particularly preferred for AGN identifications (see review by Lyu & Rieke 2022a). The galaxy dust emission template used in prospector-AGN + is based on an empirical SED template of Haro 11, which has been also tested against real observations of very highz galaxies (Lyu et al. 2016; De Rossi et al. 2018). In contrast, the other galaxy dust emission models used in other works typically do not capture some key features of highz galaxy ISM, such as low-metallicity and the possibly different dust compositions, as argued in De Rossi et al. (2018).  \nWe have normalized all models at 2 µ m. If the emission of LRDs were dominated by the dust torus at these wavelengths, considering the average redshift of our sample < z > = 6 . 5, i.e., 15 µ m observed, the slope of the torus would translate to the emission at 1 µ m rest-frame, 8 µ m observed, being ∼ 10 times fainter. The F770W-F1500W colors of LRDS are much smaller, indicating that the 1 µ m emission cannot be dominated by the torus even if the 15 µ m flux is linked to an obscured AGN.", 'C. SED FITS FOR THE FULL LRD SAMPLE IN THE SMILES/JADES FIELD': "In this Appendix, we present a detailed discussion of the SEDs of the Golden Five galaxies (shown in Figures 5 to 9) and the rest of galaxies in our LRD sample.  \nC.1. Galaxies detected up to F1800W: JADES-57356 and JADES-204851  \nFigures 5 and 6 showed the fits for JADES -57356 and JADES -204851, the 2 LRDs in our sample detected up to (at least) 18 µ m.  \nJADES -57356 is a canonical LRD, presenting a change in slope in its SED at around 2 µ m. No prominent emission lines are detected in any medium- and broad-band filter, although some excess is seen in the filters that would cover the Lyα and MgII emissions for z = 5 . 5 1 . To reproduce the characteristic, bimodal SED of the LRDs, 2 of the 4 codes use two distinct components.  \nThe SED fit with synthesizer-AGN uses 2 stellar populations (with independent attenuations), one with a young (60 Myr) and mildly unobscured (A(V) ∼ 1 mag) starburst, which also contributes to the faint emission lines, and an older (250 Myr) stellar population with much larger reddening, A(V) ∼ 4 mag, that dominates the stellar mass content (98% of the total).  \nProspector-AGN uses a similarly young (130 Myr) stellar component for the UV emission and AGN emission from a dust-obscured accretion disk for the optical. Prospector-SF , on the other hand, uses only stars to reproduce the SED up to the optical and redder wavelengths, as also obtains Prospector-AGN +. However, Prospector-SF requires an unusual, extremely gray attenuation law (n ∼ 0.4). Indeed, Barro et al. (2023) noted that a typical Calzetti law (n=0) and a single attenuation parameter would not reproduce the SED of the LRDs. In fact, as can be extracted from the data in Table 2 and the parameters written in the plots, the ratio between the far-UV (at ∼ 150 nm, FUV) and optical ( ∼ 550 nm) attenuations is ∼ 1.5 for Prospector-SF and ∼ 1.5 for Prospector-AGN +, smaller than the ∼ 2 . 6 implied by the Calzetti law. We note that the synthesizer-AGN fits assume such as law, but the combination of the independent attenuations for the old and young stellar populations results in a FUV-to-optical attenuation ratio close to the values given above ( κ FUV = 1 . 8).  \nWe also remark that synthesizer-AGN does use a Calzetti law but assumes independent extinctions for the 2 stellar populations in the host, and also for the AGN, which allows a good fit to the data for JADES -57356. ProspectorAGN + does not require a two-component fit but it allows for hybrid AGN+galaxy models, which can sometimes fit the SED with 2 distinct components, each dominating a different spectral range, but can also fit the whole SED with only one of those components. This is the case for JADES -57356 which, as noted above, exhibits a similar best-fit model to the prospector-SF model for the UV-to-NIR region.  \nOverall, all 4 models provide good fits to the UV-to-NIR SED, with some differences. The three stellar-dominated models suggest a relatively large stellar mass, M ⋆ = 10 10 . 4 M ⊙ for synthesizer-AGN . In contrast, in ProspectorAGN the emission from the obscured accretion disk (A(V) = 2 . 8 mag) dominates the SED up to λ = 1 µ m (and the AGN at even bluer wavelengths), which leads to a much smaller stellar mass for the galaxy host that is only visible in the rest-UV.  \nWhile all the models obtain similarly good UV-to-NIR fits, they differ substantially in the MIR emission, where the SED shows a clear flattening beyond 1 µ m (F770W-F1000W ∼ 0 mag) followed by an upturn around 2 µ m that steepens quickly towards the redder bands. This SED shape is well reproduced by the stellar-continuum-dominated models that combine a stellar peak around 1.6 µ m with dust emission presenting a prominent 3 µ m PAH line in the MIRI F2100W band. For synthesizer-AGN , the overall best-fitting code ( χ 2 values given in the plots), the dust emission model is quite extreme: the best template corresponds to a ULIRG with A(V) = 10 mag, 350 pc star-forming region size, 60% luminosity ratio of OB stars in hot spots compared to total luminosity, and dust density in hot spots 10 2 cm -3 . This dust emission model peaks at rest-frame ∼ 30 µ m (see Figure 16 in Appendix B), quite a blue wavelength compared to the more typical ∼ 70 -100 µ m for (U)LIRGs (e.g., Rieke et al. 2009), revealing the important role of warm/hot dust in LRDs, and also implying a relatively low emission at (sub-)mm wavelengths, where the flux of LRDs is faint (Labb'e et al. 2023b; Williams et al. 2023a).  \nThis Siebenmorgen & Krugel (2007) model points to the existence of very hot dust bathed by a very intense radiation field, coming from dust-buried OB stars (the total attenuation of stellar light in the models ranges from 2 to 4.5 mag), as indicated by the synthesizer-AGN and prospector-AGN + fits. Such an intense and compact heating source also matches well what can be expected from an AGN, which in principle could contribute to some extent to the MIR emission. Indeed, synthesizer-AGN and prospector-AGN + do show some contribution, although faint and thus uncertain, from an AGN.  \nHowever, prospector-AGN , the only code that is forced to be AGN-dominated in this spectral range, provides a worse fit to the MIRI bands because the transition from disk-dominated to torus dominated emission around λ ∼ 2 µ m leads to redder MIRI colors due to the steeper slope of the dust torus relative to the dust associated with star formation in the other models (note the bad fits to the F560W and F1000W fluxes). To some degree, prospector-SF also has problems fitting the rest-NIR region because the nebular emission partially hides the flattening of the stellar continuum and provides a worse fit to those MIRI bands.  \nIn summary, for this source the fit to the dust emission presents a smaller χ 2 when including the, most probably star-formation related, 3 µ m PAH. The relatively large stellar mass and old mass-weighted age (50-700 Myr) could support the AGN nature of the dust emission (for several reasons, e.g., the stars are relatively old and that could intuitively imply that there are not so obscured, or the evolved state of the galaxy gives time for the SMBH to grow), but the fits are significantly worse.  \nJADES -204851, whose SED and postage stamps are shown in Figure 6, is a source with spectroscopic redshift z sp = 5 . 4790 discussed previously in Matthee et al. (2023b, GOODS-S-13971 in that paper), where they reported on the presence of a broad 2200 km s -1 H α component (S/N ∼ 5 and with some artifacts in the spectrum -certainly, this is not one of the clearest BLR sources in Matthee et al. 2023b-) arising from an AGN with a 10 7 . 5 M ⊙ SMBH. This AGN would account for half of the H α emission, according to the spectroscopic analysis.  \nIn our synthesizer-AGN fits for JADES -204851 (top left panel of Figure 6), the AGN emits around 10% of the total emission at 2 µ m, also contributing to emission lines such as H α , [OIII], and MgII. The MIR emission would also have a significant contribution from star-formation powered dust emission, with a ∼ 10% possible contribution from a QSO-like emission. For this galaxy, the best fitting dust emission model corresponds to one 350 pc ULIRG star-forming region with highly embedded stars (A(V) = 72 mag), with a 40% ratio of OB stars in hot spots compared to the total luminosity, and dust density in hot spots 10 4 cm -3 . The optical-to-NIR is fitted with a combination of two stellar populations. One of them is a very young (1 Myr) unobscured (A(V) = 0 . 5 mag) starburst, which takes care of the strong emission seen in spectroscopy and detected in the NIRCam imaging and the blue continuum. The  \nother population is slightly older (10 Myr) affected by a large reddening (A(V) ∼ 2 . 5 mag) and dominating (94%) the total stellar mass (M ⋆ = 10 9 . 6 M ⊙ ).  \nThe prospector-AGN + results are very similar to those obtained with synthesizer-AGN in the rest-frame optical and near-infrared: the SED is dominated by stars and dust heated by stars. However, prospector-AGN + reproduces the UV part of the SED with stars alone, while some contribution from a (unobscured) QSO is obtained with synthesizer-AGN , mainly due to the possible presence of a MgII line at observed wavelength ∼ 2 µ m. The contribution is small, and thus uncertain, but the spectroscopy hints that there is a broad-line component (at 5 σ confidence level). This difference between both codes might also be caused by the different configurations of the stellar extinction laws. As pointed out by Kriek & Conroy (2013b), the standard Calzetti law (used by synthesizer-AGN ) typically provides poor fits at UV wavelengths for highz galaxies and thus the AGN component is likely selected from the model to fit the SED in synthesizer-AGN . Meanwhile, prospector-AGN + used the updated galaxy extinction law introduced by Kriek & Conroy (2013b), which is supposed to be more realistic. On the other hand, the independent treatment of the attenuation for old and young stars in synthesizer-AGN can overcome the problems with a single extinction parameter.  \nThe prospector-SF also provides a qualitatively good fit with a single stellar component and a modest A(V)=1.2 mag, but very gray attenuation, n = 0 . 4 (translating to κ FUV ∼ 1 . 5, same value obtained by synthesizerAGN ). Interestingly, the fit to the MIRI bands is worse due to the presence of strong emission lines in the model which would imply a larger flux in F1500W than is observed. The dust emission in this model contributes less than 10% of the total emission at 2 µ m and has virtually no impact in the best-fit SED.  \nThe prospector-AGN fit also provides a good overall fit of the NIRCam and MIRI photometry, using stars for the UV and an AGN for the optical/IR (as in the LRD characterized in Killi et al. 2023). In this galaxy, which exhibits significantly bluer MIRI colors than JADES -57356, the obscured emission from the disk (A(V)= 2 . 3 mag) dominates the optical and NIR emission up to λ ∼ 2 µ m. The emission from the torus is mostly unconstrained but it would require low IR luminosities or a very large opacity ( τ V ≳ 100). As before, the implied stellar mass of the host is the smallest of all the models, M ⋆ = 10 8 . 1 M ⊙ .  \nIn any case, the different morphology seen in F814W compared to the LW NIRCam bands, with the emission in the former being dominated by a knot located to the NW of the very concentrated (and very red) emission seen in the latter, points to star formation dominating the UV spectral range. This morphological difference between the rest-frame UV and NIR is also seen in other LRDs, for example, JADES -211388 or JADES -79803, discussed later.", 'C.2. Galaxy detected up to F1500W: JADES-219000': 'Figure 7 discussed the results for JADES -219000, the third LRD in our sample detected up to 15 µ m. This source is presented in Sun et al. (2024, in prep.), with a spectroscopic redshift of z = 6 . 8119. Similarly to JADES -57356, the best-fit with synthesizer-AGN indicates a UV-to-NIR SED dominated by stars with a relatively high mass (M ⋆ = 10 10 . 4 M ⊙ ), and dust emission revealing a very intense radiation field. The dust emission model corresponds to an intense starburst with 90% luminosity arising from OB stars embedded in a A(V) = 72 mag compact (350 pc), dense (10 4 cm -3 ) and clumpy dust cloud. The very hot dust present in this object could also be heated by an AGN, but the typical torus emission used in prospector-AGN + fails to fall within the F1800W and F2100W upper limits. This, in part, explains the larger χ 2 value compared to the other codes, jointly with discrepancies in the FUV which could be interpreted as attenuation law effects.  \nThe prospector-SF fit is also relatively good with similar best-fit values as the other galaxies, i.e., a young stellar population, with a small and very gray attenuation (A(V)=1.6 mag, n=0.33, κ FUV ∼ 2). The dust emission does not contribute significantly to the MIRI fluxes ( < 1% at λ = 2 µ m), the nebular emission dominates (but the F1000W is not well-fitted). The prospector-AGN fit also agrees well with the data. As before, the best-fit model also implies a low stellar mass for the host and a moderate attenuation for the disk (A(V)=1.6 mag), but it favors a more luminous torus with large opacity, τ V = 100. This means that the torus emission starts to contribute significantly at λ ∼ 1 µ m, but it peaks at longer wavelengths ( λ = 30 µ m) than the low opacity tori (see Figure 16 of the appendix), and therefore the model does not over-estimate the F1800W and F2100W upper limits.', 'C.3. Galaxies detected up to F1280W: JADES-211388 and JADES-79803': "Figure 8 showed the SED fits for JADES -211388, a spectroscopically confirmed LRD at z = 8 . 3846. This is one of the highest redshift LRDs and, consequently, the MIRI detections probe shorter rest-frame wavelengths than in previous galaxies.  \nThe synthesizer-AGN best-fit indicates that the SED is dominated by young (3 Myr), slightly extincted (A(V) ∼ 1 mag) stars in the blue, with some contribution from slightly older (5 Myr) and more obscured (A(V) = 2 . 5 mag) stars in the red, for a total stellar mass of M ⋆ = 10 9 . 1 M ⊙ . Remarkably, the NIR emission is dominated by a heavily obscured AGN as opposed to the star-formation powered dust emission of previous galaxies.  \nThe obscured AGN template is very similar to the best-fit model with prospector-AGN + and resembles the AGN-dominated fits in the NIR of Barro et al. (2023) or Labb'e et al. (2023b), but with significant contribution from stars up to rest-frame wavelength ∼ 2 µ m, where the AGN would contribute with 50% of the total flux), but with a lower luminosity that is still consistent with the upper limits.  \nThe prospector-SF best-fit parameters are almost identical to the previous galaxies with a young stellar population obscured with a mild but gray dust attenuation (A(V)=1.7 mag, n = 0 . 39, κ FUV ∼ 1 . 4, still failing to reproduce the FUV) and a relatively flat NIR emission dominated by the nebular continuum in the λ ∼ 1 -2 µ m range.  \nprospector-AGN also provides a good fit to the optical-to-NIR SED but it has some issues fitting the sharp Lyman break indicated by the blue NIRCam bands (we remind the reader that the attenuation is fixed to 0 mag in these models). The best-fit properties are again consistent with a young, low-mass galaxy host M ⋆ = 10 8 . 7 M ⊙ and an obscured AGN, A(V) = 3 . 4 mag. The lack of MIRI constraints at long wavelengths leads to best-fit models dominated only by the accretion disk emission even up to λ =4 µ m. However, a more significant contribution from the torus emission at λ ≳ 2 µ m would still be consistent with the upper limits, as shown in the synthesizer-AGN and prospector-AGN + fits.  \nFinally, Figure 9 shows the SED fits for JADES -79803, another spectroscopically confirmed LRD, this time at z = 5 . 4007. The SED exhibits a flattening around λ ∼ 1 µ m, probed by F560W and F770W, followed by a steep rise in F1000W and F1280W and then another flattening indicated by the upper limits in the reddest bands. The best-fit parameters with prospector-SF are similar to JADES -211388, but the stellar mass is a bit smaller, M ⋆ = 10 8 . 6 M ⊙ , with all stars being younger than 20 Myr and presenting different attenuations between 0.5 and 2 mag. The near-tomid IR emission is reproduced again by an obscured AGN, but in this case, the F1280W detection provides stronger evidence of an upturn.  \nThe prospector-AGN + also reproduces the SED with a similarly young stellar population and an obscured AGN dominating the near-to-mid IR emission.  \nThe prospector-SF best-fit continues on the same trends as before, but it provides a worse fit to the reddest MIRI bands whose steep rise can not be easily reproduced with emission from star-formation heated dust.  \nprospector-AGN provides a good overall fit, the best among all codes for this source, with the same trends in previous galaxies (low-mass host and obscured disk dominating the optical emission). In this case, the additional MIRI constraints motivate a fit with a larger contribution of the torus emission to the IR emission starting around λ ∼ 2 µ m, but consistent with the upper limits.", 'C.4. Galaxies detected up to F770W and with no MIRI detections': 'Figures 17 and 18 show the SED fits for all F770W-detected sources not presented in the main text and the source with no detection in any MIRI band, respectively. Results for the 4 fitting codes described in Section 4 are provided. The F770W sample shown in Figure 17 presents more heterogeneous SEDs than the Golden Five or F1000W samples discussed in Section 4. The global slope of the SED is more constant, with the difference between the SW and LW bands being less marked than for the other LRDs. In fact, this sample includes a very blue galaxy (187025, z sp = 6 . 9076), which only entered the sample because of the enhanced F444W flux due to an emission line. The SEDs, including MIRI upper limits (remarkably, for F1000W and F1280W), are very flat in the 1-2 µ m range, completely compatible with being dominated by stars without the need of much more contribution from other emitting components. Compared to other subsamples, the F770W sources are more extincted (nearly 3 mag on average for the stars dominating the mass) and present older mass-weighted ages (nearly 10 Myr, compared to the 3-4 Myr for galaxies detected at longer MIRI wavelengths).  \nFigure 18 shows a similar result for MIRI undetected sources, the available data (limited to NIRCam) is well fitted by stars only. The upper limits for MIRI would be consistent with hot dust emission from an AGN torus, but the current data is inconclusive.  \nConsidering the full sample of 31 sources, the median/quartiles χ 2 for synthesizer-AGN are 342 800 207 , prospectorAGN + obtains larger relative values compared to the former by a factor of 2 . 3 4 . 8 1 . 3 , as is the case for prospector-SF with 1 . 2 2 . 5 0 . 6 times larger χ 2 values, and prospector-AGN , 1 . 3 3 . 7 0 . 6 .  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 17. SED fitting results for the 7 galaxies detected at F770W and not beyond (i.e., the plot does not include any of the Golden Five or F1000W galaxies discussed in the main text). SEDs are normalized to 0.7 µ m. In gray, we show the fits to each individual galaxy, and we provide an average in brown.  \n<!-- image -->', 'REFERENCES': 'Akins, H. B., Casey, C. M., Allen, N., et al. 2023, ApJ, 956, 61, doi: 10.3847/1538-4357/acef21  \n\' Alvarez-M\'arquez, J., Crespo G\'omez, A., Colina, L., et al. 2023, A&A, 671, A105, doi: 10.1051/0004-6361/202245400 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068 Astropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12, doi: 10.3847/2041-8213/acbb08 Barro, G., P\'erez-Gonz\'alez, P. G., Cava, A., et al. 2019,  \nApJS, 243, 22, doi: 10.3847/1538-4365/ab23f2  \nBarro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2305.14418,  \ndoi: 10.48550/arXiv.2305.14418 Barrufet, L., Oesch, P. A., Weibel, A., et al. 2023, MNRAS, 522, 449, doi: 10.1093/mnras/stad947 Bisigello, L., Gandolfi, G., Grazian, A., et al. 2023, A&A, 676, A76, doi: 10.1051/0004-6361/202346219 Bogd´an, ´ A., Goulding, A. D., Natarajan, P., et al. 2023, Nature Astronomy, doi: 10.1038/s41550-023-02111-9 Bosman, S. E. I., ´ Alvarez-M´arquez, J., Colina, L., et al. 2023, arXiv e-prints, arXiv:2307.14414, doi: 10.48550/arXiv.2307.14414 Boylan-Kolchin, M. 2023, Nature Astronomy, 7, 731, doi: 10.1038/s41550-023-01937-7  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 18. SED fitting results for the 12 galaxies not detected in any MIRI band. Same format as previous figure.  \n<!-- image -->  \nBrammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786  \nBrammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2011, ApJ, 739, 24,  \ndoi: 10.1088/0004-637X/739/1/24  \nBruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000,  \ndoi: 10.1046/j.1365-8711.2003.06897.x  \nBunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467,  \ndoi: 10.48550/arXiv.2306.02467  \nBurgasser, A. J., Gerasimov, R., Bezanson, R., et al. 2023, arXiv e-prints, arXiv:2308.12107,  \ndoi: 10.48550/arXiv.2308.12107  \nByler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66  \nCalzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692  \nChabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392 Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718,  \ndoi: 10.1086/309250  \nChoi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102  \nConroy, C., & Gunn, J. E. 2010, ApJ, 712, 833,  \ndoi: 10.1088/0004-637X/712/2/833  \nConroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486  \nDe Rossi, M. E., Rieke, G. H., Shivaei, I., Bromm, V., & Lyu, J. 2018, ApJ, 869, 4, doi: 10.3847/1538-4357/aaebf8  \nDekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201,  \ndoi: 10.1093/mnras/stad1557  \nDesprez, G., Martis, N. S., Asada, Y., et al. 2023, arXiv e-prints, arXiv:2310.03063,  \ndoi: 10.48550/arXiv.2310.03063  \nMagorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353  \nMaiolino, R., Scholtz, J., Witstok, J., et al. 2023a, arXiv e-prints, arXiv:2305.12492,  \ndoi: 10.48550/arXiv.2305.12492  \n- Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023b, arXiv e-prints, arXiv:2308.01230,  \ndoi: 10.48550/arXiv.2308.01230  \n- Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023a, ApJ, 950, 67, doi: 10.3847/1538-4357/acc846\n- Matthee, J., Naidu, R. P., Brammer, G., et al. 2023b, arXiv e-prints, arXiv:2306.05448,  \ndoi: 10.48550/arXiv.2306.05448  \nNatarajan, P., Pacucci, F., Ricarte, A., et al. 2023, arXiv e-prints, arXiv:2308.02654,  \ndoi: 10.48550/arXiv.2308.02654  \nNenkova, M., Sirocky, M. M., Nikutta, R., Ivezi\'c, ˇ Z., & Elitzur, M. 2008, ApJ, 685, 160, doi: 10.1086/590483  \n- Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864, doi: 10.1093/mnras/stad2411\n- Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817  \nPacucci, F., Nguyen, B., Carniani, S., Maiolino, R., & Fan, X. 2023, arXiv e-prints, arXiv:2308.12331,  \ndoi: 10.48550/arXiv.2308.12331  \nP\'erez-Gonz\'alez, P. G., Gil de Paz, A., Zamorano, J., et al. 2003, MNRAS, 338, 508,  \ndoi: 10.1046/j.1365-8711.2003.06077.x  \nP\'erez-Gonz\'alez, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234, doi: 10.1086/523690  \n- P\'erez-Gonz\'alez, P. G., Barro, G., Annunziatella, M., et al. 2023a, ApJL, 946, L16, doi: 10.3847/2041-8213/acb3a5  \nP\'erez-Gonz\'alez, P. G., Costantin, L., Langeroodi, D., et al.  \n2023b, ApJL, 951, L1, doi: 10.3847/2041-8213/acd9d0  \n- Polletta, M., Tajer, M., Maraschi, L., et al. 2007, ApJ, 663, 81, doi: 10.1086/518113  \nRichards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470, doi: 10.1086/506525'}